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INTRODUCTION 


Observed  seismic  signals  are  superimposed  on  background  noise,  their  waveforms 
depend  on  the  characteristics  of  the  source.  They  are  also  subject  to  distortion  by  various 
propagation  effects,  and  distortion  by  the  near  source  and  recording  site  structures.  In 
nuclear  monitoring  we  need  to  discriminate  nuclear  explosions  from  earthquakes,  estimate 
yields  and  determine  other  parameters  of  interest.  Clearly,  this  requires  sophisticated 
processing  of  seismic  signals  for  reducing  the  effects  of  the  noise  and  signal  distortions  and 
to  extract  in  some  systematic  objective  fashion  the  source  parameters  we  need  to  know. 

At  the  initial  stages  of  the  VELA  Uniform  project  (1960-1970)  a  considerable 
amount  of  work  was  devoted  to  the  application  of  various  statistical  time  series  processing 
methods  to  seismic  recordings,  mostly  with  the  objective  to  reduce  the  S/N  ratio.  Such 
work  soon  led  to  disillusionment,  however,  because  many  predictions  with  regard  to  the 
improvements  in  detectability  and  the  reduction  of  noise  could  not  be  realized.  Moreover, 
the  state-of-art  of  computing  was  not  advanced  enough  for  the  large-scale  application  of 
such  methods,  computers  were  slow,  had  limited  memory  and  computer  time  was 
expensive.  Thus,  improvements  in  S/N  ratio  of  the  order  of  4-6  dB  were  hardly 
considered  to  be  worthwhile.  It  is  also  evident  now  that  many  of  the  early  assumptions 
applied  in  designing  digital  processing  schemes  were  not  correct,  most  notably  the 
assumptions  of  the  perfect  signal  coherence  and  the  stationary  background  noise.  Such 
problems  lead  to  the  virtual  abandonment  of  sophisticated  seismic  signal  processing 
methods  in  nuclear  monitoring  research,  while  the  development  continued  in  exploration 
seismology,  underwater  acoustics  and  other  fields. 
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In  addition  i<  itif  shorLComir.g  in  digital  processing  schemes  for  seismic  data,  it  was 
not  possible  to  model  the  seismograms  and  predict  the  complexities  of  the  signals 
adequately,  or  alternatively,  to  invert  the  seismograms  for  source  parameters.  For  the 
application  of  advan,;ed  signal  processing  methods,  on  the  other  hand,  we  need  realistic 
and  accurate  signal  models.  The  algorithms  for  realistic  forward  modeling  of  seismic 
waveforms  were  developed  more  recently,  and  these  require  the  speed  and  efficiency  of 
modem  computers  The  development  of  such  algorithms  brought  with  it  an  almost  totally 
deterministic  approach  to  seismic  modeling  and  interpretation  that  often  disregarded  its 
inherent  limitations  and  the  random  nature  of  data.  In  treating  seismic  data  the  only  quasi- 
statistical  methods  in  rouune  use  today  are  array  beaming  and  F-K  methods.  Source 
inversion  consists  mostly  of  direct  modeling  in  which  the  waveforms  of  raw  recordings 
are  fitted  with  deterministic  models.  Limitations  of  this  approach  are  now  quite  apparent: 
modeling  is  only  etfective  at  low  frequencies:  also  the  fact  that  most  raw  seismic  data  have 
a  very  limited  bandwidth  led  to  ambiguities  in  the  inversion  for  source  parameters.  It 
appears  that  presently  there  exists  no  comprehensive  and  objective  theory  for  broadband 
source  inversion 

It  appears  tnat  there  is  a  need  to  merge  modem  ideas  and  methods  of  signal 
processing,  estirnntion  theory,  and  seismic  modeling  with  the  purpose  of  improving  our 
capability  of  souree  identification,  location  and  detection.  Development  of  new  approaches 
to  this  problem  in  seismology  is  timely  because  of  the  advances  in  computer  technology 
and  theoretical  advancements  in  other  fields.  This  report  presents  results  and  ideas  along 
these  lines,  and  illustrates  that  such  approaches  can  be  fruitful.  Since  the  seismic  sources, 
propagation  effects,  and  signal-to-noise  ratios  can  be  best  described  as  complex  functions 
of  frequency,  and  since  the  inversion  problem  is  linear  and  thus  easier  to  solve  in  the 
frequency  domain,  most  of  our  processing  is  done  in  the  frequency  domain;  although 


2 


excursions  into  the  time  domain  are  necessary  to  interpret  deconvolution  results  and  verify 
some  of  the  basic  structures  of  signals,  such  as  spectral  factorability. 

The  report  consists  of  several  sections  of  varying  content.  The  first  major  section 
describes  the  application  and  extension  of  the  multi-channel  spectral  factorization  ideas  to 
regional  data.  First,  we  describe  a  test  that  applies  joint  deconvolution  of  regional  and 
teleseismic  data.  Following  this,  we  present  results  for  regional  Pn  and  Lg  arrivals 
recorded  at  NORESS.  We  demonstrate  that  the  concept  factorability  of  signal  spectra  into 
site  and  source  response  factors  is  applicable  for  these  phases.  Subsequently,  we  explore 
some  applications  of  such  concepts  for  event  grouping  with  respect  to  location,  source  time 
functions  and  refined  azimuth  estimation.  Central  to  such  techniques  are  the  ideas  of  inter¬ 
event  and  inter-site  transfer  functions  and  coherences. 

Following  the  application  of  spectral  factorability,  we  explore  the  possibilities  of 
new  methods  of  source  inversion,  stated  in  terms  of  extensions  and  generalizations,  of 
some  known  F-K  analysis  algonthms  (Capon,  l969,  Shumway,  1988).  Theoretical 
background  and  derivations  of  such  concepts  are  presented  in  Appendix  A  at  the  end  of  this 
report.  Specific  applications  to  teleseismic  body  and  surface  waves  are  given  in  the  main 
body  of  the  text. 

In  order  to  test  these  concepts,  fairly  extensive  simulations  were  performed  on 
synthetic  data.  Such  testing  was  also  necessary  in  order  to  debug  the  computer  codes  and 
explore  the  stability  and  inherent  limitations  of  such  methods  and  compare  them  to  other 
techniques.  In  doing  this,  we  found  that  standard  time  domain  waveform  modeling 
methods  often  do  not  fit  the  waveforms  closely  enough  to  suppon  the  conclusions  stated 
(such  as  source  finiteness  effects  and  source  depths  estimates  accurate  to  a  few  kilometers) 
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and  that  the  information  contained  in  the  seismograms  in  various  frequency  bands  is  not 
effectively  and  uniformly  utilized.  Simulations  using  synthetic  data  from  point  sources 
showed  that  the  imaging  concepts  we  derived  are  sound.  Good  source  images  were  also 
obtained  from  botli  body  and  surface  waves.  Imaging  of  simulated  faults  indicated  that  the 
method  is  suitable  for  analyzing  complex  sources. 

Following  this,  we  have  applied  the  imaging  methods  to  several  Kazakh  nuclear 
explosions  with  the  objective  to  define  the  source  depth.  The  results  showed  that  these 
sources  must  be  located  at  depths  less  than  5  km.  Applying  the  same  method  to  the  long 
period  P  waves  from  the  El  Golfo  earthquake,  gave  a  mid-crustal,  14  km  depth  that  is  not 
far  removed  from  a  previous  10  km  estimate.  Imaging  Kazakh  sources  with  surface  waves 
was  less  successful,  mostly  because  of  demonstrable  timing  errors  in  the  data  sets  we 
acquired. 

We  conclude  the  report  with  some  results  of  optimum  filtering  of  teleseismic  data  at 
EKA.  These  results  show  that  a  combination  of  site-equalization  methods  with  noise- 
adaptive  filtering  can  be  an  effective  means  for  increasing  the  signal-to-noise  ratio  of  both 
regional  and  teleseismic  data. 

Appendix  A  contains  particulars  of  the  imaging  theory  and  Appendix  B  describes  a 
straightforward  implementation  of  the  multi-channel  deconvolution  method  in  the  X- 
windows  environment.  This  implementation  opens  the  way  for  larger  scale  testing  of  this 
technique  on  teleseismic  data. 
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.1*  i^LICATION  OF  MULTI-CHANNEL  DECONVOLUTION 
AND  SPECTRAL  FACTORIZATION  CONCEPTS 
TO  REGIONAL  SEISMIC  ARRAY  DATA 

GENERAL  REMARKS 

Obtaining  source  information  from  regional  recordings  is  very  important  in  nuclear 
monitoring  because  of  the  concern  about  small,  decoupled  nuclear  explosions  possibly 
hidden  in  conventional  firing  patterns  of  quarrying  operations.  Thus  far,  little  progress  has 
been  made,  and  practically  no  work  done,  in  extracting  information  from  regional 
waveforms  by  coherent  processing,  although  spectral  expressions  of  multiple  firing 
sequences  have  been  discovered  (Willis,  J963;  Baumgardt  and  Ziegler,  1988).  The 
disadvantage  of  the  spectral  and  cepstral  methods  is  that  the  phase  information  in  the 
signals  is  discarded  and  that  the  modulation  patterns  depend  on  the  similarity  in  the 
waveforms  from  the  various  individual  shots  making  up  the  firing  sequence.  The 
advantages  of  such  methods  are  that  they  are  robust  and  do  not  otherwise  depend  on  the 
details  of  propagation.  In  this  paper  we  shall  discuss  the  general  problem  of  regional 
waveform  analysis  and  show  how  investigating  the  spectral  structure  of  the  data  may  help 
in  recovering  source  information,  in  addition  to  that  obtainable  by  other  methods  of 
analysis. 

A  general  conclusion  that  may  be  drawn  from  attempts  of  computing  synthetic 
seismograms  of  regional  arrivals  such  as  Pn,  Pg,  Sn.  and  Lg  and  comparing  them  to  real 
data  is  that  it  is  quite  impossible  to  compute  Green’s  functions  with  sufficient  precision  to 
reproduce  such  high  frequency  waveforms  in  any  detail.  The  reason  is  that  we  do  not  have 
detailed  knowledge  of  the  geological  structures  that  determine  the  waveforms,  and  even  if 
we  had  such  knowledge,  we  simply  do  not  know  how  to  compute  synthetic  seismograms 
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for  such  complicated  structures.  While  it  is  possible  to  obtain  synthetic  seismograms  for 
long  period  and  some  standard  shon  period  seismograms  that  match  observed  waveforms 
in  considerable  detail,  for  re^^ional  seisms. grams  tiie  oest  we  can  hope  for  is  to  match  the 
overall  characteristics  of  wave  envelopes  a.n.;l  reJanve  phase  amplitudes.  Details  of 
individual  waveforms  for  various  events  and  vc-cordLi.g  sites  vary  so  much,  even  for  closely 
spaced  sensors  at  small  arrays,  that  it  is  hard  to  see  any  visual  similanty  between  the  raw 
waveforms. 


Will  it  ever  be  possible  to  derive  .'Ource  rei.ned  infomtarion  from  the  aetailed 
waveforms  of  regional  seismograms  if  vse  aveep;  tiie  tact  tnat  we  shall  never  know  the 
Green's  functions?  The  answer  to  this  question,  based  on  the  work  to  be  presented  below, 
is  a  qualified  yes;  although  the  possibilities  seem  to  Pc  limited  due  to  the  extreme  variability 
of  regional  waveforms  with  source  mecha.vtsm  and  tne  positions  of  the  sources  and 
receivers. 

THEORETICAL  CONSIDERATIONS 

Before  we  review  the  results  of  our  data  analyses,  it  is  necessary  to  describe  some 
expected  spectral  characteristics  of  regional  seismic  arrivals  in  the  light  of  elementary 
seismological  theory.  Let  us  examine  the  general  spectral  structure  of  seismic  arrivals  in  an 
arbitrary,  heterogeneous  anisotropic  medium.  The  problem  is  completely  frequency- 
separable  in  the  frequency  domain  as  seen  by  transcribing  some  known  equations  (Aki  ami 
Richards,  1980;  p.  53) 

F‘io} ;  -  [  ioj (w  ;  J  ( 1 ) 


6 


into  the  frequency  domain.  St  are  the  component*  of  th  moment  tensor  and  G  are 
denvatives  of  Green's  functions  (they  will  be  called  Green's  functions  for  simplicity  in  this 
paper)  associated  with  them.  The  index  i  is  for  the  given  source  and  j  is  for  the  sensor.  In 
this  expression  the  component  index  used  by  Aki  and  Richards  (1980),  which  is  always 
vertical  in  uiis  paper,  was  dropped  and  the  summation  over  moment  tensor  components  is 
explicitly  indicated,  rather  than  implied  by  repetition  of  mdices,  in  order  to  avoid  confusion 
since  other  indices  will  be  repeated  in  some  expressions  in  this  paper.  This  structure  is 
equivalent  to  a  multi-channel  filtering  problem  in  which  we  can  consider  a  seismogram  as  a 
result  of  eitlier  of  the  following  equivalent  operations; 

a)  Filtering  the  source  inputs  for  each  component  of  a  moment  tensor 
with  filters  corresponding  to  the  Green's  functions. 

b)  Filtering  the  Green's  functions  with  the  six  independent  components 
of  the  moment  tensor. 

It  is  advantageous  sometimes  to  look  at  the  problem  using  the  second  interpretation 
since,  as  we  pointed  out  above,  we  shall  probably  never  know  the  Green's  functions  and 
any  information  we  gain  from  the  seismograms  should  depend  only  on  the  coherence 
relationships  among  the  various  seismograms.  The  relationship  (1)  is  too  general  to  be  of 
much  use,  unless  something  is  known  about  the  spectral  terms  and  factors  in  this 
expression.  First  of  all,  we  know  that  the  moment  tensor  components  generally  have  a 
short  time  response  relative  to  the  length  of  the  regional  seismograms,  since  most  seismic 
sources  we  consider  are  of  short  duration.  Moreover,  we  have  to  exploit  the  spectral 
structures  of  seismograms  to  be  discussed  below  and  may  also  invoke  the  principles  of 
reciprocity  before  we  can  derive  any  useful  information  from  regional  waveforms. 
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It  has  been  shown  that  for  shon  period  teleseismic  P  waveforms  originating  from 
sources  in  limited  regions  (such  as  nuclear  test  sites)  the  spectra  of  individual  short  period 
sensors  can  be  adequately  described  in  terms  of  the  factorability  condition  (Filson  and 
Frasier,J972 ;  Der  et  al,  1987) 


Flico)  =  s''{co)Rj{(o)  (2) 

where  R  and  5  are  the  site  and  source  factors  and  F  is  the  Fourier  transform  of  a  P  wave 
from  tl  -  source  i  and  observed  at  site  j.  The  R  are  strongly  dependent  on  the  azimuths 
(slowness  vectors)  of  the  arrivals;  thus,  this  equation  is  valid  only  for  limited  source 
regions.  We  have  been  successful  in  estimating  the  R  and  S  and  interpreting  S  for  the 
teleseismic  case  based  on  the  assumption  that  the  site  factors  can  be  reduced  by  beaming 
over  sites  and  that  the  Green's  functions  are  otherwise  simple,  and  can  be  enhanced  by  the 
same  process  (Der  et  al,  1987). 

Fortunately,  Green's  functions  for  teleseismic  short  period  P  waves  can  be 
described  in  terms  of  a  few  rays,  such  as  P,  pP  and  sP  and  possibly  a  few  multiple 
reflections  in  the  crust  near  the  source,  all  having  the  same  slownesses,  and  thus  the  source 
waveforms  can  be  interpreted  in  terms  of  a  few  phases.  This  is  not  generally  the  case  for 
regional  arrivals  which  contain  a  large  number  of  multiple  reflecting  rays  in  the  various 
crustal  layers  and  involving  various  values  of  slownesses,  even  though  a  single  regional 
phase  tends  to  be  associated  with  a  limited  range  of  slownesses.  Therefore,  even  if  we 
could  separate  the  source  and  site  terms,  they  could  not  be  interpreted  simply;  for  example, 
"pF"  would  have  different  delays  for  the  various  rays.  Alternatively,  considering  the 
regional  arrivals  as  complex  superpositions  of  high  frequency  higher  normal  mode 
wavetrains,  does  not  allow  easy  interpretation  in  terms  of  source  characteristics. 
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The  factorability  condition  given  above  suggests  a  more  useful  form  for  the 
structures  of  regional  seismograms  as  recorded  at  arrays: 

F-{(0)  =  II  Rjiw)  (3) 

P.9 


In  this  expression  we  have  factored  the  seismogram  into  spectral  factors  consisting  of  the 
moment  tensor  spectra,  factors  G  appropriate  to  the  source  region  location  and  source 
mechanism  and  finally,  the  site  factors  R  that  depend  only  on  the  type  of  wave  we  are 
analyzing  and  the  sensor  location. 

The  rest  of  this  paper  will  be  devoted  to  the  examination  of  the  properties  of  such 
factors  and  testing  of  the  validities  of  relationships  of  this  form  for  various  types  of 
regional  data.  The  form  of  relationships  given  above  includes  the  simple  factorability 
Equation  (2)  as  a  special  case.  We  must  point  out  that  despite  the  fact  that  factoring 
matrices  of  seismograms  works  for  array  recordings  of  teleseismic  signals  (Filson  and 
Frasier,  1972),  there  is  no  reason  to  believe  a  priori  that  such  factorizations  and  associated 
methodologies  will  work  for  regional  phases  which  are  very  complex  superpositions  of 
rays  or  modes,  depending  on  one's  preferred  way  of  interpreting  them.  Any  such 
factorization  must  be  tested  and  verified  by  observation  prior  to  any  attempts  to  use  them  in 
data  processing  schemes. 

We  use  the  term  "coherent  processing"  in  this  paper  for  techniques  that  exploit  the 
generalized  factorability  condition  as  defined  above.  This  is  somewhat  different  from  the 
usual  definition  coherent  processing  techniques,  such  as  beaming  and  F-K  analyses 
routinely  applied  to  regional  signals  at  seismic  arrays  (MykkelTveit  et  al,  1983;  Ingate  et  al, 
1985;  Kvaerna  and  Mykkeltveit,  1986),  which  basically  depend  on  the  assumption  of 
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complete  signal  similarity.  These  processes  and  the  regional  array  layouts  were  optimized 
for  signals  coming  from  arbitrary  directions  using  the  concept  of  averaged  signal 
coherences  vs  sensor  distances  and  are  applied  only  at  inter-sensor  distances  and 
frequencies  where  the  site  effects  can  be  neglected.  However,  if  they  satisfy  equations  of 
the  forms  (1)  and  (2),  regional  signals  across  a  regional  array  are  much  more  "coherent,”  in 
a  more  general  sense,  if  we  analyze  regional  arrivals  from  limited  source  regions.  Let  us 
consider  now  a  few  special  cases  of  Equation  (2): 

Events  With  Identical  Green's  Functions  and  Mechanisms 


This  is  the  most  commonly  studied  case,  and  numerous  simulations  of  spectral 
modulations  due  to  various  spatio-temporal  patterns  of  ripple-firing  implicitly  assume  the 
similarity  of  single  shots.  At  a  given  regional  array  far  removed  from  the  source  region, 
the  seismograms  at  the  individual  sensors  can  be  written  as 

Fi  {a))^5Uo))Hj  (CO)  (4) 

where  the  s(co)  are  the  source  functions  for  a  given  event  and  the  H(co)  are  common 
combinations  of  the  products  of  moment  tensors  and  Green's  functions  for  all  events.  A 
commonly  observed  phenomenon  for  numerous  shots  is  that  the  spectral  modulation  is  the 
same  for  all  regional  arrivals  (Baumgardi  and  Ziegler,  1988)  even  as  seen  at  several 
regional  arrays.  This  implies  a  purely  temporal  modulation,  such  as  implied  above, 
without  any  spatial  structure  for  the  source.  Scenarios  involving  combined  temporal-spatial 
shot  configurations  of  shots  can  easily  be  simulated  and  are  actually  being  applied  for 
vibration  control  in  quarrying  operations  (Anderson  and  Stump,  1989). 


Equation  (4)  above  implies  that  the  waveforms  of  events  conforming  with  this 
model  will  be  coherent  at  common  sensors  of  an  array.  Thus,  we  can  compute  the  inter¬ 
event  coherence 


C=- 


U  JJ 


(5) 


between  two  events  which  should  have  a  high  value  in  the  case  of  low  background  noise. 
Central  to  this  idea  is  the  estimation  procedure  for  the  power  spectral  compnonets  P.  These 
must  be  computed  by  some  smoothing  procedure  such  as  smoothing  the  dot  products  of  the 
signal  Fourier  transforms  denoted  as  vectors /such  that 


s 


(6) 


where  the  overbar  denotes  spectral  smoothing,  the  star  denotes  complex  conjugation  and 
the  summing  of  smoothed  products  is  performed  over  sensors.  Smoothing  can  be 
accomplished  in  various  ways:  by  actually  smoothing  the  spectra,  by  averaging  over 
shorter  time  windows,  or  some  combination  of  the  two.  In  our  case,  we  use  the 
combination  of  spectral  smoothing  and  averaging  over  array  sites.  The  equivalent  time- 
bandwidth  product  (TBWP)  of  the  result  can  be  computed  by  multiplying  together  the 
spectral  averaging  bandwidth  B,  the  total  time  length  T  of  the  data  used  for  each  array  site 
and  the  number  of  array  sites  N.  In  any  coherence  measurement  some  decision  must  be 
made  concerning  the  bandwidth  to  be  used.  High  TBWP  implies  high  stability  of  the 
coherence  results;  with  a  given  amount  of  data  this  can  be  traded  off  against  the  resolution 
in  spectral  detail  (Bendat  and  Pier  sol,  1966;  Shumway,  1988).  In  the  case  discussed  here 
we  can  assume  that  the  waveform  differences  are  associated  with  different  effective  source 
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functions,  i.e.,  firing  sequences.  These  nonnally  have  very  short  time  durations. 
Generally,  the  amount  of  smoothing  that  can  be  applied  in  coherence  calculations  depends 
on  the  time  length  of  the  impulse  response  of  the  transfer  function  between  the  two  time 
series.  As  a  rule  of  thumb,  one  can  smooth  over  a  bandwidth  of  i/T/  where  7/  is  the 

duration  of  such  an  impulse  response.  This  implies  that  heavy  spectral  smoothing  may  be 
applied  in  the  calculation  of  C  above  if  the  differences  in  waveforms  are  mostly  due  to 
different  firing  sequences.  For  further  discussions  of  the  particulars  of  coherence 
calculations  we  refer  the  reader  to  Shumway  (1988,  pages  70-73). 

Alternatively,  this  form  of  inter-relationship  implies  the  existence  of  a  transfer 
function  that  can  transform  events  one  from  the  other  with  a  filter  of  relatively  short 
impulse  response.  This  does  not  imply  waveform  similarity,  however,  although  coherence 
and  waveform  similarity  are  often  mistakenly  interchangeable. 


In  such  cases,  the  waveforms  of  individual  events  generally  will  not  be  coherent  in 
the  sense  of  Equation  (5),  and  the  spectra  of  the  superpositions  of  individual  shots  with 
different  mechanisms  wiU  not  be  modulated  in  any  easily  predictable  fashion.  A  possible 
scenario  of  such  a  situation  in  quarry  blasting  is  a  mix  of  shots,  some  of  which  may  be 
applied  to  a  vertical  rock  face  (shear  source)  and  some  buried  under  a  horizontal  surface, 
with  both  kinds  of  shots  being  close  to  the  surface  at  roughly  the  same  depth. 

Nevertheless,  as  long  as  the  six  Green’s  functions  are  common  to  a  set  of  events, 
any  seven  will  be  coherent  if  the  multiple  coherence  is  computed  for  them,  where  we  may 
test  the  possibility  of  reconstructing  the  i-th  time  series  from  the  rest  (Bendat  and  Piersol, 


1966).  We  believe  that  this  is  a  worst  case  scenario,  since  some  of  the  independent  G  will 
not  contribute  significantly  to  the  waveforms,  and  some  may  correlate  well;  thus,  only  a 
few  may  be  important.  In  such  cases,  multiple  coherences  with  lesser  numbers  of  events 
may  attain  values  significantly  different  from  zero.  No  matter  how  complex  the  source 
mechanisms  may  be,  it  will  be  possible  to  "calibrate"  parts  of  a  quarry  with  a  small  number 
of  events,  and  the  high  multiple  coherences  will  identify  all  subsequent  events  in  that 
quarry.  In  this  study  we  shall  not  present  any  multiple  coherence  results. 


This  is  analogous  to  testing  outputs  of  several  different  systems  with  independent 
random  inputs  for  coherence,  where  most  of  these  should  test  as  not  being  significantly 
different  from  zero.  This  case  covers  all  events  with  different  depths  and  at  different 
locations.  If  all  5s  are  different,  then  none  of  the  single  and  multiple  coherences  between 
that  event  and  the  events  at  a  "calibrated"  quarry  will  be  high.  This  may  be  the  case  for 
events  at  a  quarry  and  a  hidden  decoupled  nuclear  explosion  at  a  significantly  greater  depth 
or  an  earthquake  near  a  quarry.  The  lack  of  coherence  could  be  used  as  a  criterion  for  such 
events.  Given  the  fact  that  decoupled  nuclear  explosions  have  to  be  buried  deeper  than 
conventional  explosives  in  standard  quarrying  procedures,  and  that  the  effective 
mechanisms  for  the  two  types  of  events  will  no  doubt  be  different,  it  does  not  seem  to  be 
very  easy  to  mix  the  two  types  of  explosions  in  order  to  hide  testing  that  could  not  be 
revealed  by  some  kind  of  coherence  analysis. 
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‘  NALYSES  OF  DATA  FROM  QUARRY  BLASTS 


Our  data  comes  from  NORESS  recordings  of  events,  listed  in  Table  I,  identified  as 
quarry  blasts  at  the  Titania  mine  and  two  Estonian  mines  named  El  and  E4.  We  have  also 
looked  at  one  event  at  the  Blasjo  dam  site.  We  do  not  know  about  the  location  accuracy  of 
these  events  or  the  reasons,  other  than  location,  for  their  association  with  any  particular 
mine.  The  locations  of  these  events  are  shown  in  Figure  1. 

In  order  to  apply  any  of  the  ideas  above,  we  must  first  verify  the  applicability  of  the 
factorability  condition  for  any  data  set.  We  have  applied  the  multi-chaimel  deconvolution 
technique  described  by  Shumway  andDer  (1985)  and  Der  et  al  (1987)  to  two  data  sets:  the 
Pn  arrivals  from  this  group  of  Titania  mine  blasts  and  the  Lg  arrivals  from  a  set  of  mine 
blasts  in  Estonia  at  mine  E4. 
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Table  I  -  LIST  OF  “TITANIA”  AND  “ESTONIA”  EVENTS 


Event 

Date 

Origin  Time 

Latitude 

Longitude 

Magnitude 

Titania  1 

14  Feb  86 

14:13:24.9 

58.3 

6.4 

2.7 

Titania  2 

14  Feb  86* 

17:54:10.6 

58.3 

6.4 

2.3 

Titania  3 

07  Jan  86 

14:14:28.9 

58.3 

6.4 

2.2 

Titania  4 

17  Jan  86# 

14:11:01.5 

58.3 

6.4 

2.7 

E4-1 

19  Jan  88 

11:55:35.0 

59.3 

27.2 

2.5 

E4-2 

20  Jan  88 

12:23:06.0 

59.3 

27.2 

2.5 

E4-3 

21  Jan  88 

11:25:30.0 

59.3 

27.2 

2.4 

E4-4 

26  Jan  88 

09:55:37.0 

59.3 

27.2 

2.4 

E4-1 

19  Jan  88 

11:50:50.0 

59.3 

24.4 

2.5 

E4-2 

21  Jan  88 

12:53:23.0 

59.3 

24.4 

2.2 

*  Probably  removed  firom  or  not  at  the  same  part  of  quarry. 

#  Reference  event  (unmodulated). 


Starting  with  the  Pn  data  from  the  Titania  mine,  let  us  assume  initially  that  all  these 
events  had  the  same  source  nsechanism  and/or  location  and  test  the  idea  that  the  differences 
in  waveforms  are  totally  attributable  to  differences  in  firing  sequences,  i.e.,  the  factors  S 
are  common  for  all  the  events.  Inspecting  the  spectra  of  mine  blasts  from  the  Titania  mine 
in  Norway  at  NORESS,  it  can  be  observed  that  some  have  strongly  modulated  spectra, 
while  some  others  have  quite  smooth  spectra  (Figure  2)  for  all  phases,  Pn,  Pg,  Sn  and  Lg. 
This  points  to  a  common,  temporal  modulation  for  each  event,  rather  than  to  a  pattern 
imposed  by  the  spatial  configuration  of  charges  (and  possibly  also  combined  with  some 
temporal  sequence  of  explosions)  since  in  the  latter  case  we  should  see  difierent  modulation 
patterns  for  the  various  regional  phases  which  possess  notably  different  phase  velocities. 
Although  the  lack  of  modulation  does  not  rule  out  some  very  complex  spatio-temporal 
modulation  patterns  of  the  source  (although  these  must  be  quite  unlikely),  such  that  the 
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spectral  manifestations  of  it  are  suppressed,  it  is  reasonable  to  assume  that  the  spectrally 
unmodulated  quarr\'  blasts  have  simple  impulsive  source  time  functions  or  very  small 
delays  while  tlie  modulated  ones  have  been  ripple-fired.  The  Lg  phases  of  the  Estonia  mine 
blasts  are  somewhat  noisy  and  the  S/N  ratio  is  only  good  below  4  Hz  (Figures  3  and  4). 

The  decomposition  of  data  into  source  and  site  factors  is  not  totally  unique  {see  Der, 
1987),  but  by  successfully  reconstructing  the  complete  set  of  original  waveforms  from  the 
much  smaller  set  of  estimated  source  and  site  spectral  factors,  we  verified  that  the  data  has 
the  spectra'  structure  described  by  Eiquation  (2).  Figures  5  and  6  show  some  representative 
examples  of  these  reconstructions,  all  of  which  are  of  similar  quality  within  each  set. 
Typically,  the  correlation  coefficients  between  the  pairs  of  original  and  reconstmcted  traces 
were  near  or  above  0.8,  which  is  not  as  good  as  in  the  teleseismic  case,  but  still  respectable 
considering  the  high  signal  frequencies  we  utilize  here.  Note  that  the  waveforms  are  quite 
different  at  the  various  sensors  for  the  same  event  and  for  the  various  events  at  the  same 
sensors.  Joint  deconvolution  of  the  Titania  events  and  the  Blasjo  event  did  not  result  in 
acceptable  reconstructions,  showing  that  an  azimuth  difference  of  10°  is  too  large  for  using 
the  same  site  factors.  In  all  the  work  we  present  in  this  paper  we  have  used  12  sensors  of 
NORESS  roughly  evenlv  distributed  over  the  area  of  the  array. 
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Figure  4.  Representative  signal  and  noise  spectra  of  some  of  the  Lg  phases 
originating  from  Estonia  mine  events.  Only  the  Lg  phases  have  good  S/N 
ratios,  the  rest  of  them  are  burled  in  the  noise  background. 
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Figure  5  (a).  Selected  reconstructed  and  data  traces  of  Pn  arrivals  from 
the  Titania  1  event. 
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Figure  5  (b) . 


*rhfc  same  for  the  'i'ltania  event. 
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Figure  6  (a)  Selected  reconstructed  and  data  traces  of  Lg  arrivals 
from  the  E4  1  event. 
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Figure  6  (b)  The  same  for  the  E4-4  event 
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Let  us  compute  now  the  inter-event  coherences  according  to  Formulas  (5)  and  (6) 
above  for  Pn  phases  from  several  pairs  of  these  events  in  which  one  of  the  events  is  the 
relatively  unmodulated  one.  The  results  show  that  the  pairs  Titania4-1  and  Titania4-3 
(Figure  7)  show  relatively  high  average  coherences,  indicating  nearly  identical  source 
mechanism  and/or  location;  while  the  pair  Titania4-2  has  low  coherence,  indicating  that 
Event  TitaniaS  has  probably  a  different  mechanism  or  location.  Although  the  S/N  ratio  was 
the  lowest  for  Event  2,  this  coherence  still  seems  to  be  too  low  for  a  coherent  signal.  The 
TBWP  for  these  coherences  was  120. 

Let  us  hypothesize  now  that  the  unmodulated  event,  Titania4,  is  a  simple  one,  a 
single  shot.  Consequently,  the  time  domain  impulse  response  of  the  moving  average  (MA) 
filter  that  can  transform  the  waveforms  of  this  event  to  those  of  the  Titanial  must  be  a  band 
limited  representation  of  the  firing  sequence.  Therefore,  we  decided  to  design  a  Wiener 
filter  to  transform  Pn  waveforms  of  Event  Titania4  into  those  of  Event  Titanial  at  all 
sensors  of  NORESS.  Prior  to  applying  this  technique,  we  wanted  to  ensure  that  the  traces 
processed  are  as  free  as  possible  from  background  noise  and  other  distortions  of  the 
spectrum  by  applying  a  common,  minimum  phase  band  pass  filter  to  all  traces  of  both 
events  that  emphasized  the  energy  and  flattened  the  spectra  in  the  3-15  Hz  frequency  band 
where  the  signal-to-noise  ratio  was  higher. 

Computation  of  time-domain  MA  type  filters  can  be  easily  accomplished  by  several 
well  known  algorithms  using  widely  available  computer  codes  (Marple,  1987).  Time 
domain  design  is  convenient  because  we  want  to  limit  the  length  of  the  filter.  The  Pn 
arrivals  were  lined  up  for  botli  events,  windowed  and  tapered,  and  the  auto  and  cross¬ 
correlation  functions  needed  by  the  filter  design  algorithms  were  computed  by  ensemble¬ 
averaging  these  correlation  functions  over  sites.  The  ensemble-averaging  of  correlation 
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Figure  7.  Interevent  coherences  for  Pn  arrivals  from  the  Tltanla  mine. 
While  based  on  the  average  slgnal-to-nolse  ratios  most  coherences  should 
be  veil  above  .75,  the  lower  values  indicate,  for  some  events,  that  the 
signals  themselves  must  be  much  less  correlated.  The  degrees  of  freedom 
for  these  coherences  was  about  120. 


functions  over  conrmon  sites  suppresses  the  effects  of  the  site  functions  and  much  reduces 
tJiose  of  the  ambient  noise.  The  length  of  the  portion  of  the  cross-correlation  function  with 
the  highest  amplitudes  gave  a  good  indication  of  the  maximum  length  of  the  ripple-firing 
sequence,  since  the  time  length  of  the  large  amplitudes  in  this  function  cannot  exceed  this 
length.  For  the  cases  discussed,  this  upper  limit  is  roughly  of  the  order  of  0.5  seconds. 
The  compound  correlation  functions  constitute  the  inputs  to  the  Levinson  algorithm  used 
for  the  computation  of  the  cross-equalizing  Wiener  filter  consisting  of  31  weights. 

Subsequently,  the  filter  was  applied  to  the  simple  event  to  derive  waveforms  of  the 
complex  events  at  all  sites.  A  set  of  traces  showing  the  representative  improvement 
between  the  two  events  at  common  sensors  are  shown  in  Figure  8.  Note  that  the  original 
waveforms  are  much  more  different  for  the  two  events  and  how  the  waveforms  became 
more  similar.  Nevertheless,  the  increase  in  the  time  domain  correlation  coefficients  was 
only  modest,  from  0.59  to  0.75  for  the  first  pair,  and  0.49  to  0.57  for  the  second.  This 
indicates  that  either  we  need  a  longer  filter  or  that  the  event  pairs  are  too  far  from  each 
other,  in  either  case,  we  cannot  attribute  the  waveform  difference  entirely  to  the  source  time 
functions.  Tiie  same  method  did  not  work  at  all  for  the  pair  Titania4-2. 

We  have  also  applied  the  cross-equalization  for  the  whole  wavetrains,  consisting  of 
all  the  regional  arrivals  of  Events  Titania4-1.  The  processing  resulted  in  a  spectacular 
increase  of  the  correlation  coefficients  between  the  two  events  from  .08  to  .65  over  the 
whole  set  of  sensors,  again  confirming  that  the  waveforms  of  these  events  are  closely 
related. 
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Figure  8  (a)  Wiener  filtering  results  for  transforming  the  Pn  phases  of 
event  Tltanla  4  Into  those  of  event  Tltanlal. 


28 


Titanla4  -♦  TitaniaS 


B5 


TItaniaS 


Titania4 


Titania4  TilaniaS 


Trtania3 


I _ I 

2.0  SEC 


Figure  8  (b)  Viener  filtering  results  for  transforming  the  Pn  phases  of 
event  Tltanla  4  into  those  of  event  Titanla  3. 


29 


Having  established  that  such  techniques  may  be  suitable  for  identifying  groups  of 
,  Jilts  with  nearly  identical  mechanisms  and/or  locations,  we  shall  examine  the  nature  of 
:  elative  inter-site  transfer  functions  for  the  same  event.  We  know  that  such  transfer 
.  ;s  exist,  and  that  they  are  consistent  from  event  to  event.  Otherwise  we  would  not 
ii  able  to  jointly  deconvolve  and  reconstruct  all  these  events.  It  would  be  desirable 
.sfcr  functions  would  be  simple,  similarly  to  the  inter-event  transfer  functions 
kior  the  coherent  event  pair,  because  then  relative  displacement  over  distances  of  the  order 
of  the  NORESS  diameter  would  not  destroy  the  inter-event  coherence  (applying  the 
principles  of  reciprocity).  In  that  case,  the  lack  of  coherence  among  events  in  a  limited 
source  region  observed  at  an  array  would  be  mostly  due  to  differences  in  source 
mechanisms. 

U-ifortunately,  this  is  not  the  case.  We  have  tried  to  design  short,  31  point,  time 
liiicrs  for  equalizing  Pn  waveforms  at  sensors  D1  and  D5  for  the  suite  of  four 
w .  :i.i ,,  bui  they  essentially  failed  to  produce  even  remotely  similar  waveforms.  Thus,  the 
inter- site  transfer  functions  are  generally  complex,  with  very  long  time  domain  impulse 
responses.  Applying  the  reciprocity  argument,  this  implies  that  a  comparably  small 
displacement  in  the  location  of  a  source  would  totally  destroy  the  waveform  similarity, 
ssuming  that  the  source  region  and  the  NORESS  site  are  similar  in  the  degree  of 
jiOgeneity  in  structure. 


Site  Equalization  Processing 

A  major  problem  in  locating  regional  events  is  that  the  first  arrivals  are  often  small 
and  buried  in  noise  and  that  most  phases  are  emergent,  without  clearly  defined  arrival 
times.  Azimuths  estimated  from  F-K  analyses  are  often  not  very  accurate,  since  the  site 
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effects  tend  to  destroy  the  plane  wave  character  of  the  signals,  thus  broadening  the  main 
lobes  of  array  response  patterns.  There  are  conflicting  demands  of  keeping  the  array 
apertures  small  for  ensuring  signal  similarity  and  increasing  the  directional  resolution  of  the 
arrays.  The  fact  that  the  site  transfer  functions  are  consistent,  albeit  not  simple,  still  leaves 
open  the  possibility  of  utilizing  this  internal  consistency  for  refining  the  estimation  of 
relative  azimuths  among  closely  spaced  events.  If  the  forms  of  the  inter-site  transfer 
functions  do  not  change  much  with  small  azimuth  changes,  except  for  some  small  time 
shifts,  then  such  a  property  could  be  exploited.  To  test  this  idea  we  have  band-pass  filtered 
the  four  Titania  events  to  emphasize  the  band  with  the  maximum  signal-to-noise  ratio, 
between  3-10  Hz.  Subsequently,  we  have  taken  the  frequency  domain  site  factors  derived 
from  the  multi-channel  deconvolution  calculations  and  multiplied  the  Fourier  transform  of 
trace  D1  with  the  site  factor  of  D5  and  vice  versa.  This  is  basically  an  inter-correlation 
technique  applied  to  the  broadband  site  (rather  than  event)  equalization.  The  resulting 
spectra  were  then  multiplied  in  the  frequency  domain  and  transformed  back  into  the  time 
domain,  as  described  by  the  formula 


<I>it)=F 


(9) 


to  view  the  cross-correlation  function.  In  this  expression  the  d’s  are  the  Fourier  transforms 
of  the  traces  and  the  R’s  are  the  site  transfer  functions.  The  results  for  three  of  the  events 
are  shown  in  Figures  9  and  10;  Event  2  was  too  noisy  and  was  left  out.  In  inspecting  these 
figures,  we  must  point  out  that  the  results  of  cross-equalization  do  not  have  to  look  like 
normal  seismograms,  with  clearly  defined  phases.  In  multiplying  with  the  complex  site 
factors,  we  essentially  perform  a  circular  convolution  with  considerable  wraparound  effect. 
Figure  8  demonstrates  that  the  site-equalized  trace  waveforms  indeed  became  quite  similar. 
Any  small  shifts  in  the  peaks  do  not  seem  to  be  associated  with  noise  or  mismatch  since  all 
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Figure 

sensors 


Titania  1  corr.  coeff.  *  0.82 


Titania  3  corr.  coeff. »  0.56 


Titania  4  Corr.  coeff. »  0.78 


9.  Cross  correlations  of  the  site  equalized  traces  of  Pn  between 
Dl  and  D5  for  selected  Titania  events. 
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Site  equalized  Pn  traces  for  selected  Titania  events. 
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cross-correlation  functions  for  this  set  of  events  have  a  well  defined,  single  maximum  with 
little  ambigutity  in  the  time  of  the  maximum. 

The  same  process  was  repeated  with  the  recordings  of  Lg  phases  from  quarry  blasts 
in  Estonia  as  recorded  at  NORESS.  Applying  the  process  to  four  E4  mine  events,  we 
obtained  good  trace  equalizations,  with  all  correlation  peaks  lined  up  at  zero  lag.  Using 
four  events  from  E4  and  two  from  El  mines  in  a  joint  site-equalization  process,  we  have 
found  that  the  efficiency  of  the  inter-site  equalization  drastically  deteriorated  as  indicated  by 
the  decrease  in  the  correlation  coefficients  between  equalized  traces  (Figure  11).  The 
correlation  peaks  are  still  generally  at  the  same  times,  but  one  of  the  peaks  (for  Event  E2)  is 
shifted  in  time.  The  decrease  in  site-equalization  efficiency  may  be  attributed  to  the  fact  that 
mines  E4  and  El  are  too  far  from  each  other  to  have  identical  site  functions  at  NORESS. 
The  waveform  similarity  is  not  even  visible  for  this  set  (Figure  12). 

Given  the  fact  that  the  waveforms  become  quite  similar  after  applying  some  cross¬ 
equalization  treatment  to  the  outputs  of  a  pair  of  sensors  located  at  the  extremities  of  a  small 
array  for  events  in  some  limited  source  areas  (not  the  two  Estonia  mines  together),  it  may 
be  not  justifiable  to  limit  the  sizes  of  regional  arrays  to  a  few  kilometers.  It  appears  that 
regional  signals  from  limited  source  regions  are  actually  “coherent,”  in  the  sense  of  a 
broader  definition  of  this  term,  over  much  larger  distances  between  sensors.  Clearly,  the 
maximum  array  diameters  dictated  by  the  average  waveform  correlation  properties 
(Mykkeltveit  et  al,  1983;  Der  et  al,  1988)  do  not  apply  here.  It  could  be  fruitful  to  test  and 
develop  similar  techniques  using  data  from  small  arrays  larger  than  NORESS. 
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Figure  11.  Cross  correlations  of  the  site-equalized  traces  of  Lg  bcLwcen 
sensors  Dl  and  D5  for  selected  events  at  the  Estonia  mine  E4. 
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JOINT  DECONVOLUTION  OR  TELESEISMIC  AND  REGIONAL  DATA 


Given  the  complex  nature  of  regional  arrivals,  it  would  be  very  fortuitous  if,  similar 
to  teleseismic  data,  we  would  be  able  to  see  simple  P,  pP  or  sP  arrivals  after 
deconvolution.  Unlike  teleseismic  arrivals  which  can  be  described  in  terms  of  a  direct  P 
wave  and  several  distinct  surface  reflections,  regional  arrivals  are  made  up  of  numerous 
generalized  rays  with  various  slownesses.  It  is  indeed  not  possible  to  see  anything  similar 
to  P,  pP  and  sP  sequences  in  multi-channel  deconvolutions  of  regional  arrivals.  A  logical 
modification  of  the  teleseismic  multi-channel  deconvolution  procedure  for  incorporating 
regional  data  could  consist  of  using  only  the  teleseismic  data  for  constructing  the  source 
function  estimate,  thus  defining  a  regional  site  ftinction  which,  after  dividing  it  out,  would 
yield  the  usual  simple  sequence.  In  in  order  to  test  this  idea,  we  have  identified  two 
Kazakh  events  for  which  EKA,  GBA  and  KAAO  data  were  available.  In  this  set  the  station 
KAAO  contributed  the  regional  data  (Pn),  and  the  two  UK  arrays  the  teleseismic  data. 
Joint  deconvolution  of  all  these  (using  the  scheme  in  Figure  13)  showed,  using  the  events 
listed  in  Table  n,  that  the  factorization  worked;  i.e.,  the  reconstructions  were  of  good 
quality  at  both  at  the  teleseismic  UK  array  and  the  regional  KAAO  sensors  (Figure  14). 
This  seems  to  indicate  that  schemes  like  this  could  work  on  larger  data  sets.  Nevertheless, 
drawing  major  conclusions  from  data  sets  this  small  would  be  premature.  It  was  also 
difficult  to  find  data  where  we  had  both  teleseismic  and  regional  data  of  good  quality,  i.e., 
where  neither  clipping  occured  or  the  teleseisms  were  not  swamped  by  background  noise. 
Widescale  testing  and  possible  applications  for  such  schemes  may  have  to  wait  until 
systems  with  very  large  dynamic  ranges  and  32  bit  digitizers  will  become  common. 
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TABLE  n  -  KAZAKH  EVENTS  USED  FOR  JOINT  DECONVOLUTION 
OF  REGIONAL  AND  TELESEISMIC  DATA 


Event 

Origin  Time 

Lat  (N) 

Lon  (E) 

1978241 

04:26:75.9 

49.82 

78.14 

1978162 

02:56:57.6 

49.90 

78.80 
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Figure  13.  An  iterating  scheme  to  jointly  deconvolve  teleselsmlc  and  regional 
P  waves.  Only  the  teleselsmlc  data  are  used  for  estimation  of  source  factors 
but  all  traces  are  utilized  for  deriving  the  site  responses. 


Multichannel  Deconvolution  Dla 


Figure  14.  Reconstructions  of  some  of  the  data  at  ERA  and  GBA 
(teleselsmlc) ,  and  at  KAAO  (Pn— regional)  from  the  spectral  factors 
estimated  by  the  scheme  In  Figure  13. 


BROADBAND  IMAGING  OF  SEISMIC  SOURCES 
AND  SOURCE  INVERSION 


GENERAL  REMARKS 

This  part  of  the  report  presents  a  new,  general,  statistics-based,  broadband 
frequency  domain  approach  to  source  inversion  that  avoids  much  of  the  subjective  bias  that 
is  prevalent  in  time  domain  direct  modeling  and  source  inversion  methods.  This  method  is 
also  very  suitable  for  automation  of  the  source  inversion  and  seismic  discrimination 
process. 

MOTIVATION  FOR  FREQUENCY  DOMAIN  INVERSION 

A  popular  and  widely  used  method  for  the  derivation  of  source  and  path  parameters 
from  teleseismic  recordings  is  the  direct  modeling  of  the  observed  seismic  waveforms.  The 
objective  of  the  procedure  is  to  obtain  some  fits  between  the  observed  waveforms  and  the 
synthetic  seismograms  optimized  in  some  sense.  At  this  point  it  is  assumed  that  the 
parameters  used  in  computing  the  synthetic  seismograms  are  the  best  estimates  for  those  of 
the  source  and  path.  One  can  make  some  interesting  observations  reading  such  papers. 
The  most  obvious  is  that  the  synthetic  and  real  seismograms  do  not  match  perfectly  and  that 
the  differences  are  usually  in  the  finer,  high  frequency  details  of  the  waveforms. 
Moreover,  the  synthetic  seismograms  usually  contain  less  high  frequency  than  the  data 
even  if  we  allow  for  the  presence  of  background  noise.  The  inverted  source  time  functions 
generally  contain  fine  details  that  depend  on  the  same  high  frequencies  that  were 
mismatched  in  the  waveform  modeling  process.  Numerous  publications,  using  long  period 
waveforms  as  data  claim  to  have  determined  source  depth  to  an  accuracy  of  a  few 
kilometers  and  have  seen  indications  of  source  finiteness.  Such  claims  in  many  cases  seem 
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almost  too  good  to  be  true,  considering  the  wavelengths  of  the  seismic  waves  at  the 
frequencies  where  the  waveforms  seem  to  fit.  Most  studies  use  unmodified  data  for  fitting 
that  were  obtained  from  systems  with  narrowband  characteristics,  such  as  those  of  the 
WWSSN  network.  The  qualifiers  "generally"  and  "most"  are  used  above  to  signify  that 
these  observations  are  not  valid  for  all  time  domain  work  of  this  type  and  the  quality  of  fits 
and  the  claims  made  also  vary  considerably.  Neverthelesss,  because  of  the  widespread  use 
of  such  techniques  such  apparent  discrepancies  should  be  of  concern. 

While  it  is  easy  to  make  these  observations  visually,  we  need  some  measures  of  the 
misfits  between  synthetic  seismograms  and  data.  Let  us  quantify  first  the  general  standards 
for  time  domain  fitting  of  long  period  data.  To  obtain  some  quantitative  measure  of  the 
generally  accepted  standards  of  waveform  fits  we  have  taken  22  synthetic  and  data 
waveforms  pairs  randomly  selected  from  the  literature,  enlarged  and  digitized  them  on  a 
digitizing  table.  The  differences  between  the  data  and  the  synthetic  traces  are  clearly  visible 
in  the  original  publications,  and  are  not  products  of  the  digitizing,  i.e.,  these  differences  are 
much  larger  than  any  plausible  digitizing  errors. 

In  order  to  obtain  a  measure  of  goodness-of-fit  in  the  time  domain  we  computed  the 
correlation  coefficients  between  data  and  synthetic  waveforms 


This  is  also  the  most  common  measure  for  judging  the  quality  of  fit  in  time  domain 
waveform  inversion  studies  (if  any  quantitative  measures  of  waveform  similarity  are 
applied  at  all).  For  our  selected  sample  a  histogram  of  correlation  coefficients,  shown  in 
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Figure  15,  illustrates  that  the  average  time  domain  fits  we  have  chosen  have  correlation 
coefficients  averaging  near  0.9  with  a  few  over  0.95.  Several  of  the  waveform  pairs  are 
also  shown  in  Figure  16.  These  figures  should  convince  the  reader  that  we  have  not 
purposely  chosen  poor  waveform  fits  to  make  our  point,  but  that  these  are  better  than  the 
average  fits  in  the  literature.  The  synthetic  seismograms  also  contain  less  high  frequency 
than  the  data  in  accordance  with  visual  observations  that  anyone  can  make  by  reading  such 
studies  (Figure  17).  Moreover,  the  spectra  are  sharply  peaked,  giving  the  data  a 
narrowband  character.  We  would  like  to  establish  some  relationship  between  the  quality  of 
time  domain  fits  and  the  resolution  in  the  spatio-temporal  configuration  of  sources  derived 
fi'om  the  data. 

Let  us  evaluate  now  the  data  vs  synthetic  fits  in  the  frequency  domain,  using  FFT’s 
of  the  data  and  synthetic  traces,  and  using  an  alternative  similarity  measure  for  the  same 
data  set.  This  measure  is  an  ensemble-averaged  coherence 


di  ico) 5* {(0) 


C(cu)  = 


Vs  d.  (co)d.*((u)  s-  (CO)  s*  ico) 


summing  (ensemble  averaging)  over  waveform  pairs  i  and  applying  spectral  smoothing  as 
denoted  by  the  overbar.  By  plotting  this  as  a  function  of  frequency  (Figure  18),  the 
problem  with  time  domain  waveform  fitting  with  narrowband  data  becomes  obvious. 
While  the  fit  near  the  peaks  of  the  spectra  are  good  (but  by  no  means  perfect)  with  a 
coherence  a  little  over  0.9,  this  fit  rapidly  decreases  in  quality  with  increasing  fiequency. 
Although  we  are  confident  that  our  digitization  is  good  at  least  to  0.3  Hz,  the  coherence  at 
this  frequency  and  below  it  is  essentially  indistinguishable  from  zero  for  this  data  set.  Note 
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Figure  15.  Histogram  of  time  domain  correlation  coefficients  between 
synthetic  seismograms  and  data  randomly  selected  from  the  literature. 
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Figure  17.  Comparison  of  averaged  spectra  of  selected  data  and 
s^mthetic  seismograms. 


Figure  18.  Ensemble-averaged  coherences  of  our  sample  of  data-synthetic 
pairs  defining  the  "average  fit." 
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that  this  test  could  not  even  be  performed  with  data  from  any  single  study  since  practically 
none  of  them  use  as  many  as  22  stations.  We  shall  designate  this  coherence  plot  as 
rep’^esenting  the  quality  of  a  "standard  fit"  in  time  domain  modeling  of  long  period  data. 
Clearly  this  standard  will  have  a  profound  effect  on  the  resolution  of  the  time-space  details 
of  sources. 

With  regard  to  the  source  time  function  it  should  be  obvious,  when  looking  at  the 
inverse  source  problem  formulated  in  the  frequency  domain,  that  it  is  simply  not  possible  to 
say  anything  about  the  higher  frequency  details  of  the  source  function  at  all  if  the 
corresponding  frequency  components  are  not  matched  in  the  data  as  indicated  by  Figure  16. 
This  is  because  the  problem  becomes  completely  frequency-separable  in  the  frequency 
domain  as  seen  by  transcribing  the  well  known  equations  from  Aki  and  Richards  (1980) 

into  the  frequency  domain 

u„{x,0))  =  M^^  i(o).G^^io))+N„  {(0) 

The  problem  in  the  frequency  domain  also  becomes  much  simpler  and  easier  to  solve 
(Olson  and  Anderson,  1988)  and  instead  of  taking  the  time  domain  trace  amplitudes  as 
data,  we  could  fit  the  amplitudes  and  phases  of  the  various  frequency  components  at  the 
various  recording  stations.  It  also  follows  by  simply  inspecting  the  formulation  above  that 
descriptions  of  fine  details  found  in  the  space-time  structure  of  seismic  sources  are  not 
credible  unless  the  quality  of  time  domain  fits  is  improved  at  higher  frequencies.  This  is  a 
very  obvious  and  trivial  point  that  is  frequently  overlooked. 
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The  points  made  above  give  a  basic  reason  for  our  search  for  alternative  methods 

and  some  quantitative  tests  for  source  inversion,  and  we  return  to  these  arguments  after  the 

next  section.  We  shall  show  that  the  average  quality  of  time  domain  fits  as  described  by  the 

coherences  in  Figure  18  are  barely,  or  not  at  all,  sufficient  for  deriving  source  depth  and 
« 

other  source  parameters  to  any  useful  accuracy.  In  fact,  most  work  of  this  type  seems  to 
operate  on  the  border  of  non-resolvability  and  marginal  significance.  In  the  next  section  we 
describe  a  new,  promising  alternative  approach  to  source  inversion.  We  do  not  claim  that 
our  approach  is  the  best  possible  one;  other  methods,  just  as  good  or  better,  could  probably 
be  found. 

METHODS  FOR  SEISMIC  IMAGING 


Looking  at  the  mathematical  problem  of  using  a  set  of  seismic  signals  from  a 
number  of  seismic  stations  and  trying  to  estimate  the  energy  coming  from  a  given  location 
of  the  source  region,  it  can  be  shown  to  be  identical  in  mathematical  form  to  F-K  analysis 
where  the  vectors  of  beam  steer  operators  (phasors)/are  replaced  by  some  more  general 
operators  (Green's  functions).  Instead  of  having  components  with  absolute  value  unity,  as 
in  the  array  beam-steer  case,  now  we  may  have  to  normalize  the  "beam"  results  by  dividing 

with//*.  Nevertheless,  we  have  as  before  a  number  of  F-K  algorithms  available  to  us  in 
this  application.  The  first  one  is  the  generalized  beam,  robust  and  plagued  by  side  lobes: 


the  second  is  the  Capon  (1969)  estimator,  a  high  resolution  operator  with  poor  stability 

c=[/-s-vi 
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And  finally,  the  Shumway  estimator  a  good  compromise  between  the  first  two; 

fSL 


ff 


where  S  is  the  spectral  matrix  of  the  input  data  (typically  an  outer  product  of  the  Fourier 
transforms  with  a  small  proportion  of  unit  matrix  added,  see  Appendix  A).  Using  various 
appropriate  /  vectors,  we  can  now  examine  the  source  radiation  in  space  (horizontal 
coordinates  and  depth).  We  shall  discuss  the  properties  of  these  statistics  in  some  detail  in 
Appendix  A. 

Ideas  similar  to  these  were  advanced  by  Goldstein  and  Archuleta  (1987)  in 
seismology  and  Baggeroer  et  al  (1988)  in  underwater  acoustics  using  the  MUSIC  algorithm 
and  the  Capon  algorithm  respectively  (Capon,  1969).  The  third  algorithm  we  prefer  does 
not  have  some  of  the  undesirable  properties  of  Capon's  method  in  this  application.  The 
statistic  has  an  approximate  noncentral  F  distribution  at  individual  frequencies  under  the 
null  hypothesis  (no  signal)  with  2  and  2(N-1)  degrees  of  freedom  for  input  waveforms. 
If  we  stack  the  results  for  the  same  source  coordinates  obtained  at  various  frequencies,  the 
probability  density  function  of  the  result  is  approximately  normal  with  mean  Mm  and 
variance  with  M  the  number  of  frequencies  used  in  the  stack,  and  m  and  o’ being  the 
mean  and  the  standard  deviation  of  the  individual  F  distributions. 
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APPLICATIONS  OF  THE  FORMALISMS  TO  VARIOUS  SEISMOLOGICAL 
PROBLEMS 


Let  us  assume  now  that  we  have  observed  seismic  body  waves  from  a  seismic 
source  at  numerous  stations  at  various  distances  and  azimuths  from  the  source.  In  this 
case,  we  know  the  wavenumbers  at  each  station  for  each  particular  wave  type  may 
construct  an  image  of  the  source.  Assuming  that  we  know  the  source  mechanism  of  an 
elementary  source  (slip  on  an  element  of  a  fault  plane  at  a  given  depth)  and  that  the  source 
region  is  horizontally  homogeneous  (flat  layered),  we  may  write  the  components  of  vector 
/  as 


f^  (m)  =  C,.  (0))  exp  (j  (0  s-  x  ) 

where  the  first  factor  Ci(co)  is  a  frequency  domain  radiation  pattern  appropriate  to  the  crust 
and  to  the  azimuth-slowness  for  P  waves  seen  at  station  i,  i.e.,  the  sequence  of  P  waves 
following  the  various  types  of  rays  possible  for  a  source  of  that  type  at  the  given  depth;  the 
exponential  factor  corresponds  to  horizontal  translation  with  respect  to  the  source  position 
vector  X,  and  the  unit  wavenumber  vector  appropriate  to  a  given  station  i  is  denoted  by  k- 
In  laterally  heterogeneous  structures  the  translation  must  be  described  by  more  complicated 
functions,  but  the  general  principles  are  the  same. 

Alternatively,  we  may  wish  to  analyze  surface  waves  observed  at  a  set  of  stations  in 
a  similar  way.  Then  the  components  of  the  / vector  become 

/;  {(0)=P.i0))exp[j  (ok^  jf/c,  (G))] 
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In  this  case,  the  waves  are  dispersive  and  the  phase  velocity  depends  on  frequency,  and 
Pi( (o)  is  the  source  radiation  pattern  in  the  direcdon  of  the  station  /.  The  exponential  again 
describes  horizontal  translation  in  laterally  homogeneous  models. 

It  is  easy  to  see  that  this  approach  can  be  generalized  to  source  region  structures  of 
any  complexity  and  any  combinations  of  observations.  For  laterally  heterogeneous 
structures  no  simple  exponential  factors  could  be  used  for  effecting  horizontal  translation, 
of  course,  but  the  essential  mathematical  forms  shown  above  would  still  be  valid. 

The  formalism  also  provides  an  objective  way,  as  opposed  to  visual  fits  in  the  time 
domain,  to  compare,  over  a  chosen  frequency  band,  synthetic  seismograms  to  data.  In  this 
application  the  5  can  be  an  outer  product  of  observation  vectors  (FFT’s  of  data)  and  the/s 
could  be  the  vectors  of  FFT's  of  synthetic  seismograms  using  various  parameterizations. 
The  parameterizations  that  maximize  any  of  the  statistics  defined  above  wiU  be  the  best  fits. 
This  way  many  of  the  problems  with  time  domain  fitting  discussed  below  can  be  avoided. 

We  present  three  figures  to  illustrate  that  these  concepts  work  and  that  the  statistics 
indeed  have  maxima  at  coordinates  (and  other  parameters  as  well)  which  coincide  with  the 
actual  positions  of  sources.  In  Figure  19  we  show  a  source  image,  in  the  form  of  depth 
slices  of  a  point  source  located  at  6  km  depth  at  the  center  of  the  2  dimensional  grid.  This 
image  was  constructed  using  Shumway's  formula.  In  Figure  20  we  show  some  synthetic 
Rayleigh  waves  which  were  used  for  constructing  the  image  in  Figure  21.  The 
corresponding  image,  constructed  with  Capon's  formula,  is  shown  in  Figure  21. 
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0.05H2<f  <0.333  Hz 


depth  m  12  km 


depth  m  15  km 


Figure  19.  Inages  of  a  double-couple  point  source  at  various  assumed 
source  depths  with  5Z  %fhlte  noise  constructed  from  synthetic  P  wave 
seismograms  at  10  teleselsmlc  stations  well  distributed  In  azimuth 
and  slowness.  This  S/N  ratio  corresponds  to  a  fit  with  a  correlation 
coefficient  near  0.99.  Shumway's  algorithm  was  used  In  this  example. 
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Figure  20.  Image  of  a  double-couple  point  source  with  5%  white  noise 
constructed  from  synthetic  Rayleigh  wave  seismograms  at  8  teleseismic 
stations  well  distributed  in  azimuth.  Capon's  algorithm  was  used  in 
this  example. 
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Figure  21.  Synthetic  surface  waves  used  for  constructing  the  image 
in  Figure  20. 
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SIMULATION  RESULTS 


The  Capon  and  Shumway  algorithms  have  been  implemented  in  FORTRAN  and 
tested  for  correcmess.  These  computer  programs,  besides  their  uses  in  data  analyses,  are 
also  suitable  for  investigating  the  sensitivity  of  these  methods  to  perturbations  in  the  data 
and  for  the  testing  of  the  resolution  obtainable  from  such  methods  for  various  source 
parameters,  most  notably  in  the  spatial  details  of  the  sources.  In  the  following  we  present 
the  results  of  some  simulations. 


Having  defined  the  statistics  providing  source  image,  let  us  examine  how  misfitting 
waveforms  in  the  fashion  described  in  Figure  18  will  distort  a  source  image,  and  affect  the 
results.  For  simplicity  we  shall  assume  the  source  crust  to  be  a  halfspace  so  that  the 
generalized  beamsteer  vector  becomes 


where  the  subscript  I  refers  to  the  particular  station,  and  the  subscripts  j=l,2^  refer 
to  the  phases  P,  pP  and  sP,  respectively. 

Where  ifc  is  a  unit  vector  appropriate  to  a  station  i,  x  is  a  position  vector  in  the 
source  region,  the  c/i  are  phase  velocities  at  a  given  station,  the  tij  are  the  travel  times  of  the 
P,  pP  and  sP  phases.  The  source  image  is  a  mapping  of  F  over  the  source  region  in  depth 
and  horizontal  coordinates.  Depth  dependence  enters  into  this  formulation  by  the  delay 
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nines  of  the  depth  phases  and  the  exponential  factor  will  be  determined  by  the  horizontal 
coordinates  of  the  source.  We  thus  assume,  for  the  sake  of  convenience,  that  the  Green's 
functions  are  perfectly  known,  which  is  generally  not  the  case.  The  common  problem  of 
poorly  defined  Green's  functions,  which  makes  the  inverse  problem  even  less  tractable,  is 
not  addressed  here. 

The  statistic,  F,  becomes  infinite  for  the  case  of  perfect  fit,  i.e.,  in  the  noiseless 
case.  Having  developed  a  statistically  tractable  representation  of  the  source,  we  apply  this 
algorithm  now  to  a  typical  inversion  problem,  one  with  ten  azimuthally  well  distributed 
stations  (30  degrees  apart)  and  having  three  associated  P  wave  slownesses  (0.05, 0.06  and 
0.08  sec/km).  The  simulated  data  consisted  of  teleseismic  P  waves  from  a  double  couple 
source  (strike=245°,  dip=30°  and  rake=120°)  situated  at  6  km  depth.  In  Figure  19  we 
have  shown  images  of  this  source,  located  at  the  center  of  the  search  grid,  at  different 
depths.  In  this  figure  we  have  increased  the  trace  term  in  Shumway's  formula  by  5%  in 
order  to  keep  the  statistic  finite  and  used  the  frequencies  between  .05  and  .333  Hz  with 
equal  weight. 

To  investigate  how  this  clear  image  degrades  as  we  perturb  the  signals,  we  decrease 
the  coherent  part  of  the  inputs  by  some  factor  and  add  some  uncorrelated  random 
contributions  to  them  such  that  the  average  power  at  each  frequency  remains  roughly  the 
same;  but  the  average  coherence  at  that  frequency  achieves  the  desired  value.  Using  these 
corrupted  inputs  in  the  imaging  process,  we  can  assess  whether  the  misfitted  synthetic 
waveforms  are  still  defining  the  same  source  image  as  the  data,  which  is  assumed  to  be 
correct.  Alternatively,  we  may  consider  these  images  as  the  true  source  entirely  consistent 
with  the  perturbed  inputs  and  the  original  clear  images  as  alternatives  to  them  that  fit  the 
data  within  some  tolerances. 
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In  the  first  simulation  we  have  added  about  10%  uncorrelated  noise  (in  amplitude 
ratio)  evenly  throughout  the  frequency  band  quoted  above.  This  corresponds  to  about  a 
0.99  correlation  coefficient  between  the  original  and  perturbed  waveforms  having  flat, 
broadbpd  characteristics  in  the  band  indicated.  These  waveforms  still  give  a  somewhat 
rugged-looking,  but  clear,  well-defined  image  of  our  source  (Figure  22).  We  must  note, 
however,  that  correlation  coefficients  between  the  synthetic  and  original  waveforms 
associated  with  this  problem  are  far  better  than  those  commonly  seen  in  waveform  smdies. 

Next  the  standard  for  the  "average  fit",  defined  by  the  coherences  in  Figure  18,  is 
applied  to  the  problem.  Again,  we  add  a  proportion  of  random  uncorrelated  noise  to 
reproduce  the  coherences  for  "average  fit"  and  construct  an  image.  We  would  expect  that 
for  acceptable  misfits  (perturbations)  the  main  features  of  the  solution  would  remain  intact. 
This  expectation  is  not  met,  however,  since  it  can  be  seen  in  Figure  23  that  the  image 
disappeared  totally.  Repeating  the  same  procedure  in  a  narrower  frequency  band  of  0.05- 
0.15  Hz,  in  order  to  utilize  only  the  frequency  components  which  are  correlated  with  the 
highest  coherences  in  Figure  18  (near  0.9)  in  Figure  24,  we  still  have  some  very  poorly 
defined  peaks.  The  maximum  is  at  the  wrong  depth  and  is  off  center  and  there  are  maxima 
of  roughly  the  same  size  at  greater  depth  near  the  edges  of  the  grid. 

As  a  final  test,  we  decided  to  reverse  the  trend  to  progressively  pooer  fits  and 
increase  the  S/N  (coherent  and  random  noise  ratio)  in  the  narrower,  0.05-0.15  Hz, 
fi^uency  band  to  conespond  to  a  conelation  coefficient  of  about  0.97.  With  this  value  we 
are  starting  to  see  the  image  emerging  again,  but  the  maxima  at  3  and  6  km  of  depth  are  the 
same  to  four  decimals  (we  marked  them  both  in  the  associated  Figure  25).  On  tfte  positive 
side,  the  maxima  at  greater  depth  are  considerably  smaller;  thus  it  can  be  claimed  that 
the  depth  was  resolved  to  an  accuracy  of  ~6  km.  Nevertheless,  to  achieve  this,  it  was 
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Figure  22.  Image  %rlth  lOZ  random  noise  added,  otherwise  the  synthetic 
data  used  were  the  same  as  in  Figure  19. 
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Image  of  the  same  point  source  by  adding  noise  such  that  the 
corresond  to  the  "average  fit"  in  the  0.05-0.333  band. 
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Figure  24.  Image  of  the  same  point  source  by  adding  noise  such  that  the 
coherences  correspond  to  the  "average  fit"  in  the  0.05-0.15  band. 
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Figure  25.  Image  of  the  same  point  source  by  adding  noise  such  that  the 
coherences  correspond  to  the  time  domain  correlation  coefficient  of  .97 
In  the  0.05-0.15  band. 
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necessary  to  increase  the  degree  of  correlation  to  a  level  not  commonly  attained  in  data- 
synthetic  comparisons.  It  is  disturbing  to  find  that  the  imaging  procedure  breaks  down  at 
precisely  the  correlation  level  deemed  "acceptable"  by  most  time  domain  studies. 

The  above  findings  do  not  imply,  of  course,  that  the  conclusions  drawn  from  all 
direct  waveform  modeling  studies  are  invalid.  Many  studies  do  not  make  excessive  claims 
with  regard  to  fine  details  of  source  time  functions,  source  depth,  and  faulting  history. 
Many  are  no  more  than  point  source  solutions  where  some  surface  reflections  were  utilized 
to  refine  the  source  parameters.  The  qualities  of  waveform  fits  and  the  number  of  stations 
utilized  also  vary  considerably,  although  generally  most  fits  are  no  better,  and  many  are 
much  worse,  than  the  ones  we  picked  for  analysis.  Nevertheless,  the  instability  of  the 
results  with  respect  to  small  perturbations  in  the  fits  of  the  various  frequency  components 
of  the  signals  indicate  that  many  published  conclusions  in  such  studies  must  be  marginal  at 
best  and  probably  could  not  be  substantiated  even  by  some  simple  elementary  statistical 
tests.  It  appears  that  although  comparisons  of  synthetic  and  data  waveforms  may  appear  to 
be  impressive,  they  do  not  seem  to  be  sufficient  to  ensure  stable  or  even  valid  results  and 
that  good  visual  waveform  fits  and  high  synthetic -data  correlation  coefficients  of  narrow- 
band  raw  seismograms  mean  much  less  than  it  is  generally  assumed. 

Simulating  a  Fault 

The  examples  shown  above  illustrate  that  the  method  works,  i.e.,  images  have 
maxima  where  the  actual  sources  are  located.  Instead  of  imaging  point  sources,  it  is  more 
interesting,  firom  a  practical  point  of  view,  to  examine  how  cases  with  multiple  sources  can 
be  handled.  In  Figure  26  we  show  an  image  of  a  fault,  made  up  of  a  row  of  point  sources 
diagonally  arranged  in  the  x-y  plane  at  the  same  depth.  This  image  was  formed  by  using 
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Fault  trace 


100  km 


Figure  26.  A  fuzzy  image,  indicative  of  the  resolution  that  may  be  achieved  with 
good  quality  real  data  for  large  events,  of  a  simulated  fault  obtained  using 
Shumway's  algorithm.  The  fault  runs  at  a  45“  angle  with  the  sides  starting  from 
the  center  of  the  plot. 
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Shumway's  method  and  using  the  /  vectors  appropriate  to  a  single  point  source.  In 
imaging  we  have  inceased  the  trace  of  the  matrix  in  the  equation  for  Shumway’s  statistic  by 
a  factor  of  2.  This  improves  the  image  by  reducing  a  large  maximum  generated  by  the  low 
frequency  components  of  the  signal  that  otherwise  would  dominate  the  image.  This  image 
is  intermediate  between  a  generalized  beam  and  the  one  that  we  would  get  from  the  exact 
Shumway  formulation.  Simulations  like  these  are  ideal  for  assessing  the  resolution  of  the 
method  and  its  dependence  on  the  number  of  stations,  their  azimuthal  and  distance 
distribution,  and  the  frequency  range  utilized.  The  poor  resolution  in  Figure  26  is  the  result 
of  the  fact  that  only  eight  simulated  stations  were  used.  Although  it  may  not  be  obvious, 
such  limitations  are  equally  valid  for  other  types  of  source  inversion  studies. 


Another  interesting  problem  is  that  of  the  separation  of  explosion  and  double  couple 
contributions  in  a  set  of  seismic  surface  waves.  An  obvious  application  is  yield  estimation 
from  surface  waves  when  tectonic  strain  release  is  present.  We  want  to  reduce  the  effects 
of  the  double  couple  component  and  estimate  the  energy  in  the  compressional  source 
contribution.  For  this  application  the  Capon's  algorithm  seems  to  be  the  most  suitable, 
since  it  is  essentially  a  formulation  for  the  energy  in  an  output  of  a  constrained  optimum 
filter  designed  to  reduce  the  noise  (double  couple  component)  and  pass  the  signal 
(compressional  component)  undistorted.  The  original  formulation  of  Capon's  F-K 
algorithm  has  the  inverse  noise  spectral  matrix  N  instead  of  that  of  data  S.  In  most  imaging 
applications  (Baggeroer  et  al,  1988)  of  this  the  total  data  spectral  matrix  containing 
contributions  from  both  the  desired  signal  and  the  undesirable  noise  (strain  release  in  our 
case)  is  substituted.  This  may  have  serious  consequences,  however,  in  the  case  of 
mismatch  (Cox,  1973).  In  our  simulations  of  extracting  the  compressional  contribution  we 
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shall  use  the  original  formulation  of  Capon  (1969)  and  will  use  the  spectral  matrix 
constructed  only  from  a  priori  estimates  of  the  average  strain  release  sources  assumed  to  be 
thrust  events  (typical  for  Kazakh).  The  resulting  algorithm  will  be  more  stable  and  much 
less  sensitive  to  mismatch.  It  must  be  possible,  however,  to  obtain  some  estimate  of  the 
average  strain  release  spectral  matrix,  although  this  estimate  need  not  be  very  accurate. 

Table  HI  shows  the  results  of  such  a  simulation  where  a  constant  compressional 
source  superimposed  on  variable  levels  of  a  strain  release  source  (thrust  fault)  was 
processed  by  Capon's  algorithm.  The  "noise"  spectral  matrix  was  constructed  of  a 
superposition  of  thrust  fault  contributions  (that  also  included  the  simulated  fault  used  in  the 
input)  the  fault  orientations  of  which  were  varied  over  a  range  of  five  degrees  to  account  for 
the  uncertainties  and  variations  in  the  strain  release  components.  The  table  shows  that  the 
power  in  the  output  (signal  estimate)  did  not  vary  much,  although  the  ratio  of  the 
compressional  and  "strain  release"  components  did.  Only  at  extremely  large  ratios  can  we 
see  a  breakdown  of  the  process.  This  indicates  that  such  a  processor  could  work,  despite 
the  uncertainties  and  variations  in  the  strain  release  components.  In  contrast  to  the 
presently  used  painstaking  procedures  of  fitting  radiation  patterns  to  the  data,  estimating 
and  subtracting  the  strain  release  components,  this  procedure  can  be  automated,  performed 
instantaneously  after  receiving  the  surface  waves  from  the  event. 
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TABLE  in  -  OUTPUT  OF  CAPON’S  ALGORITHM 
FOR  SOURCE  SEPARATION  EXPERIMENTS 


Relative  Earthquake  Component  Output  Amplitude 


0. 

1.0 

0.17 

1.0 

0.333 

0.99 

0.5 

0.99 

1. 

0.99 

5. 

1.1 

10. 

1.4 

The  only  problem  that  needs  to  be  solved  is  to  obtain  in  practice  good  prior 
estimates  of  the  compressional  and  strain  release  spectral  matrices.  These  can  be  obtained 
empirically,  either  by  finding  events  with  low  strain  release  (small  Love  waves)  or 
manipulating  several  events  with  similar  strain  release  mechanisms,  but  different  degrees  of 
strain  release.  Such  a  procedure  involves  translating  the  two  events  to  the  same  location  by 
phase  shifting  and  differencing  their  spectra  to  obtain  estimates  of  a  pure  compressional 
source.  This  process  would  totally  eliminate  the  need  for  complex  “path  corrections.”  It  is 
possible,  however,  that  the  method  may  not  work  over  larger  areas  (exceeding  70  km  in 
diameter)  due  to  waveform  distortions  caused  by  scattering  of  Rayleigh  waves  (Rivers  and 
von  Seggern,  1979).  This  would  not  detract  from  the  usefulness  of  such  an  approach, 
since  other  methods  such  as  that  of  Stevens  (1986)  would  be  also  affected  by  the  same 
problem. 
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DATA  ANALYSES 

Imaging  of  Selected  Nuclear  Explosions  at  Kazakh. 

In  order  to  test  our  scheme  of  source  imaging  we  have  selected  two  Kazakh  events 
for  which  we  had  deconvolutions  at  three  UK  arrays  (Der  ex  al,  1987).  Taking  the 
deconvolved  traces  at  these  three  arrays,  we  have  used  the  model  of  an  explosion  source  in 
a  halfspace  to  design  the /  vectors  in  the  procedure.  Since  we  have  only  three  waveforms, 
complete  imaging  in  the  horizontal  plane  was  not  practical  since  it  is  not  possible  to  make  a 
good  image  with  only  three  azimuths.  Therefore,  we  have  tried  to  match  the  depth  only. 

We  must  point  out  that  while  we  are  trying  to  match  an  ideal  model,  the  actual 
deconvolutions  may  not  contain  a  true  pP  in  the  sense  of  an  elastic  model.  Moreover,  even 
if  the  simple  elastic  model  were  applicable,  the  bandwidth  limitations  of  the  data  may  not 
allow  us  to  resolve  the  true  depth.  These  caveats  notwithstanding,  if  we  plot  Shumway’s 
statistics  as  functions  of  source  depth  they  indicate  shallow  source  depths  (Figures  27  and 
28).  Thus,  the  deconvolved  seismograms  best  resemble  an  elastic  compressional  source  at 
shallow  depths;  this  is  more  than  we  can  deduce  from  location  calculations  at  teleseismic 
distances. 


It  appears  that  our  imaging  procedure  may  be  a  good  equivalent  or  substitute  for  the 
discrimination  and  depth  estimation  schemes  devised  by  Pearce  (1980)  who  also  tried  to 
match  the  various  depth  phases  with  amplitude  limits  obtained  from  measurements  made  on 
seismic  records.  Our  approach  differs  from  his  in  that  we  can  use  the  whole  seismograms 
automatically  without  manual  measurements;  we  also  utilize  the  polarities  of  phases,  but 
without  trying  to  determine  them  subjectively,  and  that  we  can,  in  addition,  incorporate  any 
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Figure  27.  Shuinway's  statistic  as  a  function  of  depth  computed  from 
deconvolved  P  wave  seismograms  of  the  December  25,  1975  Kazakh  event. 


Figure  28.  Shumway's  statistic  as  a  function  of  depth  computed  from 
deconvolved  P  wave  seismograms  of  the  June  29,  1977  Kazakh  event. 
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available  knowledge  of  crustal  structures,  no  matter  how  complex,  into  the  scheme  by 
using  them  to  compute  the /vectors.  The  matching  procedure  will  work  in  spite  of  any 
complexities  in  waveforms  caused  by  non-impulsive  source  time  functions.  Our  new 

standard  method  of  discrimination  will  require  that  we  line  up  the  P  waves  on  the  first  onset 

« 

(the  only  manual  step),  invoke  the  best  crustal  response  of  the  source  region  and  match 
these  against  sets  offs  computed  for  explosion  sources  at  various  depths  in  that  structure. 
In  more  general  crustal  structures,  complexities  not  explainable  to  simple  depth  phases  will 
also  be  present,  and  these  can  also  be  matched.  Good  matches  (high  F  values)  will 
designate  explosions,  poor  matches,  various  crustal  earthquakes  (because  the  mismatches 
in  the  amplitudes  ana  polarities  of  the  various  arrivals,  mostly  of  the  depth  phases). 
Deconvolution  is  not  needed  for  the  success  of  the  scheme,  because  common  frequency 
factors  cancel  in  the  process;  on  the  other  hand,  corrections  for  indqual  path  attenuation  or 
instrument  responses  will  be  required.  Reduction  of  site  effects  by  array  averaging  is 
desirable  prior  to  this  process,  but  does  not  seem  to  be  absolutely  necessary  if  many 
stations  are  being  used.  Crustal  structures  have  been  the  subject  of  studies  over  many 
decades  and  have  been  mapped  all  over  the  world.  A  data  base  of  such  structures  could  be 
a  good  accessory  to  intelligent  monitoring  systems. 


^ 11  !!  I  SI  IwF' x? 


In  order  to  apply  our  algorithm  to  real  data,  we  tried  to  find  some  simple  event  that 
would  be  easy  to  model  and  match.  To  avoid  computing  and  correcting  for  site  effects,  we 
decided  to  use  long  period  data,  where  such  effects  can  be  neglected,  at  least  in  the  first 
approximation.  The  El  Golfo  earthquake  was  reported  to  be  a  simple  point  source,  a  strike 
slip  event  found  to  be  at  10  km  depth  by  Ebel  et  al  (1978).  The  source  time  function  was 
estimated  to  be  a  simple  triangle  of  4  sec  duration  in  the  same  study,  and  the  t*  along  the 
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paths  to  the  receiving  stations  recording  on  WWSSN  15-100  LP  instruments  was  assumed 
to  be  1.3  sec.  We  adopted  fault  orientation  and  the  crustal  model  given  in  the  same  paper 
and  used  Shumway's  formula  in  an  attempt  to  verify  the  depth  estimate  using  the  P  wave 
data. 


The  Green's  functions  (f  vectors)  were  calculated  by  using  an  adaptation  of  the 
formulation  of  body  wave  radiation  from  a  point  source  given  by  Douglas  et  al  (1971) 
generalized  to  a  flat-layered  source  crust  model.  The  computer  program  in  question,  named 
PSV,  was  written  by  the  first  author  of  this  report  during  the  early  seventies.  The  synthetic 
seismograms  generated  (Figure  29)  are  very  close  to  the  ones  given  in  the  paper  by  Ebel  et 
al  (1978).  We  attribute  some  small,  occasional  differences  to  the  fact  that  while  our 
program  automatically  includes  all  rays  in  the  source  crust,  the  program  used  by  Ebel  et  al 
needs  the  specification  of  all  rays  desired  (Langston,  1976)\  thus,  some  less  important  rays 
may  have  been  omitted  in  the  calculations  of  Ebel  et  al.  In  any  case,  the  synthetic 
waveforais  we  have  generated  are  matching  the  data  just  as  well  as  those  of  Ebel  et  al.  We 
must  note,  however,  that  we  do  not  consider  these  fits  very  good;  they  are  worse  than 
those  used  for  computing  the  coherences  in  Figure  18. 

The  long  period  P  waveforms  to  be  processed  were  obtained  by  digitizing  the 
enlarged  seismograms  from  the  paper  on  a  high  resolution  digitizing  tablet.  Compared  to 
the  sizes  of  the  mismatches  between  synthetics  and  data,  the  digitizing  errors  and  trace 
distortions  must  be  very  small  and  inconsequential,  and  should  not  affect  the  conclusions 
of  this  test.  In  order  to  eliminate  the  absolute  times  as  a  factor  in  the  matching  procedure  (if 
we  had  known  the  absolute  times,  properly  corrected  for  mantle  heterogeneity,  we  could 
have  obtained  the  depth  from  those  alone),  we  have  lined  up  both  the  synthetics  and  the 
data  on  the  first  P  arrival.  We  assumed,  first,  that  the  relative  amplitudes  are  appropriate  to 
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El  Golfo  earthquake. 


S^..thct:lc  seismograms  appropriate  to  10km  depth  for  the 
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the  radiation  patterns  in  both  the  synthetics  and  the  data  in  the  paper.  We  have  also  made 
computations  normalizing  both  the  data  and  synthetics  to  unit  maximum  amplitudes.  In  the 
latter  case  we  used  all  the  data  traces  with  equal  weight,  which  is  justifiable  since  the 
background  noise  is  small  and  the  traces  with  small  amplitude  carry  just  as  much 
information  about  the  source  as  the  ones  with  larger  amplitudes.  In  fact,  some  of  these, 
recorded  at  the  stations  OTT  and  WES,  were  crucial  in  the  original  study  in  constraining  the 
fault  planes. 

We  have  made  two  test  runs  while  varying  the  source  depth,  fitting  the  spectral 
matrix  structures  in  the  .02  to  .  1 5  Hz  band,  using  all  frequencies  with  equal  weight.  In  the 
first  we  matched  the  data  with  the  source  crustal  responses  to  sources  at  different  depths;  in 
the  second  we  matched  the  synthetic  seismograms  for  10  km  depth  with  source  crustal 
responses  at  other  depths.  As  one  would  expect,  the  latter  resulted  in  a  maximum  of  the 
statistic  at  10  km  depth.  The  data,  on  the  other  hand,  do  not  fit  a  source  at  10  km  depth 
when  we  apply  the  same  procedure  to  them.  The  broadband  fits,  with  and  without 
normalization,  seem  to  require  a  depth  near  14  km  (Figures  30  and  31).  This  is  not  much 
of  a  difference,  but  it  is  considerably  outside  of  the  2  km  accuracy  limits  claimed  by  Ebel  at 
al  (1978)  for  their  solution.  Despite  the  fact  that  the  original  El  Golfo  earthquake  study  also 
utilized  some  additional  data,  a  limited  amount  of  SH  and  surface  wave  data,  we  do  not 
believe  that  the  inclusion  of  these  would  have  changed  our  results  or  conclusions.  We 
prefer  our  own  solution,  because  it  is  based  on  a  broadband  fit,  and  because  it  is  based 
solely  on  the  structure  of  the  data  spectral  matrix,  without  tradeoffs  against  the  source  time 
functions  and  associated  parameters  that  may  occur  in  global  waveform  modeling,  but 
which  are  factored  out  in  our  procedure.  It  is  gratifying  to  find,  however,  that  we  were 
able  to  obtain  a  solution  in  the  fi-equency  domain  for  a  mid-crustal  event  similar  to  those 
found  by  other,  more  detailed  studies.  Nevertheless,  the  task  of  testing  the  statistical 
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Figure  30.  Shumway's  statistic  for  the  synthetic  seismograms  (top) 
and  the  long  period  P  wave  data  for  the  El  Golfo  earthquake  (bottom) 
as  functions  of  depth.  Trace  amplitudes  were  not  normalized. 
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Figure  31.  Shumway's  statistic  for  the  synthetic  seismograms  (top)  and  the 
long  period  P  wave  data  for  the  El  Golfo  earthquake  (bottom)  as  functions 
of  depth.  All  trace  amplitudes  were  normalized  to  unity. 


significance  of  this  difference  still  remains  to  be  done.  The  sizes  of  differences  between  the 
depth  estimates  obtained  by  visual  fitting  of  time  domain  waveforms  and  our  broadband 
procedure  are  not  surprisingly  large  in  view  of  the  simulations  shown  earlier  in  this  report. 

Analyses  of  Surface  Waves  From  Kazakh  Nuclear  Explosions 

During  this  phase  of  the  work  we  have  collected  surface  wave  recordings  of 
Kazakh  nuclear  explosions  at  common  stations  using  the  data  bases  at  the  Center  for 
Seismic  Studies.  We  have  collected  surface  wave  data  for  nine  Kazakh  explosions  which 
were  recorded  at  five  common  stations  ANTO,  GRFO,  MAT,  TATO  and  KONO.  The  first 
test  was  to  verify  that  the  events  gave- maxima  near  to  their  actual  locations.  This  was 
found  to  be  the  case  for  several  of  them.  In  Figure  32  we  show  the  maps  of  Capon's 
statistics  obtained  by  using  the  event  with  the  smallest  tectonic  release,  lowest  Love  wave 
energy,  as  a  master  event  defining  the /vector.  This  seems  to  confirm  the  findings  of  von 
Seggern  (1972),  who  was  able  to  locate  events  relative  to  each  other  with  surface  waves 
with  about  a  10  km  accuracy.  However,  some  other  events  did  not  locate  well  at  all.  Upon 
examining  the  data,  we  have  found  that  the  absolute  times  given  in  the  data  base  were 
grossly  in  error  for  at  least  two  of  the  traces.  Further  inquiries  revealed  that  people  who 
had  used  the  data  before  (Richard  Baumstark  for  instance)  also  experienced  problems  with 
the  absolute  times  in  some  of  the  CSS  data  bases.  While  gross  errors  are  easy  to  notice, 
smaller  errors  may  go  unnoticed  and  can  make  any  results  questionable.  Therefore,  we 
decided  that  we  need  to  verify  absolute  times  for  all  records  we  shall  use  in  the  future 
before  we  attempt  anything  like  the  separate  estimation  of  compressional  and  deviatoric 
source  energy. 
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Figure  32.  Location  results  obtained  by  imaging  surface  waves  for  two 
Kazakh  explosions  (October  12,  1980  on  top  and  September  14,  1980  on 
the  bottom)  by  the  Capon  algorithm.  The  geographical  coordinates  of  the 
centers  of  the  print  plot  grids  are  50N-79E  and  the  events  are  plotted 
relative  to  that.  Filled  circles  denote  the  locations,  obtained  from 
body  waves,  as  listed  by  unclassified  files  at  the  Center  for  Seismic 
Studies,  triangles  show  the  maxima  of  the  Capon's  statistics.  The  grids 
are  50  km  on  each  side,  and  the  tops  are  facing  North. 
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ADAPTIVE  FILTERING  STUDIES 

GENERAL  REMARKS 

In  applying  discrimination  and  yield  estimation  procedures  to  measurements  made 
on  various  telesesismic  and  regional  wave  arrivals  emitted  by  the  sources,  it  is  important  to 
extend  the  ranges  of  detectability  and  measurability.  Commonly,  large  events  will  excite 
body  waves,  Lg,  and  long  period  surface  waves  that  can  be  seen  teleseismically.  Smaller 
events  will  not  have  long  period  surface  waves  that  can  be  seen  above  the  noise  at 
teleseismic  distances,  but  will  still  have  detectable  P  waves.  Events  that  are  even  smaller,  t, 
can  be  detected  only  by  regional  stations  and  arrays.  With  decreasing  magnitude,  the 
quantity  and  quality  of  measurements  used  as  discriminants  and  yield  estimators  will  also 
decrease.  Therefore,  we  need  techniques  that  extend  the  ranges  of  detectability  of  various 
weak  arrivals  and  that  enable  us  to  extract  broadband  information  from  the  data.  A  family 
of  such  techniques  can  be  loosely  characterized  as  adaptive  filtering  methods.  A  common 
trait  of  many  of  these  is  that  they  were  designed  to  suppress  noise  while  passing  coherent 
signals  without  distortion  (Booker  and  Ong,  1971;  Shen,  1979). 

Adaptive  filtering  of  array  data  was  performed  previously  by  Shen  (1979),  who 
applied  maximum-likelihood  adaptive  filters  to  short  period  teleseismic  P  waves  at  small 
arrays.  He  found  that  the  noise  reduction  exceeded  that  of  the  conventional  beamforming. 
The  filters  were  designed  on  the  basis  of  perfect  similarity  of  seismic  signals.  This 
assumption  is  not  correct  for  most  arrays,  as  we  have  pointed  out  above.  Therefore,  we 
have  investigated  the  effect  of  increasing  the  signal  similarity  by  factoring  out  the  site 
responses  prior  to  applying  the  technique.  The  factoring  out  of  site  responses  will,  of 
course,  change  the  noise  coherence  structure  also,  but  we  rely  on  the  adaptation  of  the  filter 
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to  adjust  to  this  changed  structure.  After  removing  the  site  responses,  we  can  assume  with 
more  justification  that  the  signals  are  perfectly  coherent.  This  will  reduce  the  "beam-loss" 
to  much  lower  levels. 

TELESEISMIC  P  WAVES  FROM  KAZAKH  RECORDED  AT  EKA 

For  testing  this  approach,  we  have  used  recordings  of  a  set  of  Kazakh  (Shagan 
River)  P  waves  recorded  at  EKA  and  estimated  the  site  effects  by  applying  our 
deconvolution  procedure  to  the  seismograms  listed  in  Table  IV.  Subsequently,  we  factored 
out  the  site  effects  and  applied  Shen's  algorithm  to  the  result.  In  order  to  test  the  effect  of 
the  algorithm  on  the  noise  alone,  we  have  applied  the  noise-adaptive  filter  to  the  noise 
preceding  the  signals.  We  performed  this  test  on  noise  samples  with  and  without  the  site 

TABLE  IV  .  LIST  OF  KAZAKH  EVENTS 
RECORDED  AT  AWRE  EKA  ARRAY 


Event 

Origin  Time 

Lat  (N) 

Lon  (E) 

1978241 

04:26:75.9 

49.82 

78.14 

1978162 

02:56:57.6 

49.90 

78.80 

1977334 

04:06:57.4 

49.96 

78.89 

1976342 

04:56:57.4 

49.96 

78.85 

effect  correction.  Figures  33  to  35  show  that  noise  adaptation  reduces  the  noise  amplitudes 
relative  to  the  conventional  beam  throughout  the  frequency  range  of  0-4  Hz,  but  the  noise 
reduction  is  most  effective  below  1  Hz.  The  removal  of  site  effects  did  not  impair  the 
effectiveness  of  noise-adaptation.  The  filter  adapts  quickly  (as  evidenced  by  a  few  large 
oscillations  in  the  output  after  the  noise  amplitude  falls  off). 
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Figure  33.  Comparisons  of  noise  adaptation  without  (bottom)  and  with 
(top)  the  site  response  removed.  Comparing  the  shapes  and  levels  of 
spectra,  it  can  be  seen  that  the  overall  noise  is  significantly  reduced 
by  the  adaptive  beaming  method  relative  to  the  beams,  especially  below 
2  Hz.  The  adaptive  beams  quickly  adapt  to  the  noise,  and  after  a  few 
oscillations  the  overall  noise  levels  are  much  reduced  relative  to  the 
beam  outputs. 
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Figure  34.  Teleselsmlc  P  waves  processed  by  beaming  and  adaptive  filtering 
with  and  without  removal  of  site  frequency  responses. 
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Figure  35.  Spectra  representative  to  the  signal  beam  loss  (top)  and  the  S/N 
gain  obtained  for  the  two  events  by  using  the  four  modes  of  processing 
described  In  the  text. 


Applying  the  process  to  windows  containing  signals  shows  that  the  noise  was 
reduced  to  a  level  below  visibility  in  this  high  S/N  sample  by  the  adaptive  process,  while 
the  noise  is  clearly  visible  on  the  beams  (top  traces  in  Figure  34).  Obviously,  this  example 
has  no  practical  use,  since  the  signal  is  much  larger  than  the  noise  in  the  case  shown.  It 
must  be  pointed  out,  however,  that  the  relative  noise  reduction  is  the  same  regardless  of  the 
relative  signal  size  (the  site-equalized  process  will  be  unaffected  by  the  presence  of  large 
coherent  signals).  Thus,  a  weaker  signal  would  have  stood  out  much  better  on  the  trace 
processed  by  the  same  procedure  than  on  the  beam  and  the  result  would  be  much  better 
suited  for  any  further  analyses  of  the  signal. 

It  is  also  interesting  to  look  at  these  noise  reduction  results  in  the  tiequency  domain. 
In  Figure  34  we  show  the  differences  in  the  output  on  signal  windows  when  the  adaptive 
beam  algorithm  is  applied  with  and  without  site  effect  correction.  The  figure  shows  that 
without  the  site  effect  correction  the  signal  output  is  considerably  less  with  adaptive 
beaming  than  with  conventional  beaming  (a).  This  is  not  surprising,  since  the  adaptive 
beam  process  will  reduce  the  signal  output  if  the  signals  are  not  matched.  After  site- 
equalization  is  performed,  this  signal  loss  is  much  less  severe  (b).  Since  the  noise 
reduction  of  adaptive  beaming  remains  more  effective,  even  after  applying  site-equalization, 
than  the  conventional  beam,  the  resulting  overall  S/N  ratio  improvement  will  be  better  for 
the  adaptive  beam  combined  with  site-equalization  (Figure  35).  Thus,  this  combined 
process  appears  to  be  superior  to  the  conventional  beaming  process.  The  gain  is  achieved 
at  a  considerable  expense  in  processing  complexity,  but  it  may  be  justified  in  special  cases 
when  we  want  to  learn  more  about  smaller  events  at  teleseismic  distances.  Such 
procedures  may  be  part  of  a  toolbox  for  more  detailed  investigations. 
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WESTERN  NORWAY  EARTHQUAKES  RECORDED  AT  NORESS 


Similarly  to  teleseismic  signals,  we  have  done  a  limited  amount  of  work  testing  the 
suitability  of  adaptive  beaming  and  site-equalization  methods  for  enhancing  regional  Pn 
arrivals.  Pn  arrivals  are  generally  emergent  at  regional  distances  and  are  much  smaller  in 
amplitude  than  Lg  or  Sn.  Enhancing  Pn  arrivals  may  be  quite  useful  in  discrimination 
studies.  The  methodology  followed  is  identical  to  that  described  for  teleseismic  P 
elsewhere  in  this  report.  The  method  was  tested  using  a  set  of  Pn  arrivals  from  western 
Norway  earthquakes.  A  list  of  these  events  is  given  in  Table  V.  The  site  effects  were 
estimated  from  applying  the  factorization  procedure  using  the  multi-channel  deconvolution 
program  (Der  et  al,  1987).  Applying  the  adaptive  beaming  program  to  the  noise  alone 
(Figure  36)  shows  that  the  site  response  removal  somewhat  degrades  the  relative 
performance  improvement  due  to  the  adaptive  beaming  process.  On  the  other  hand,  the 
decrease  in  signal  loss  through  the  process  appears  to  compensate  approximately  for  this 
loss.  The  net  result  is  that  the  best  strategy  seems  to  be  to  apply  adaptive  (not 
conventional)  beaming  either  with  or  without  site  corrections.  We  must  add  that  this  test 
was  done  before  we  realized  the  extreme  sensitivity  of  site  factors  to  azimuth,  and  the 
western  Norway  earthquakes  may  be  too  scattered  in  azimuth  to  get  stable  and  effective  site 
correction  factors.  In  any  case,  the  adaptive  beaming  process  appears  to  be  better  than 
conventional  beaming  for  regional  signals. 
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TABLE  V  -  LIST  OF  WESTERN  NORWAY  EVENTS 
RECORDED  AT  NORESS 


Event 

Origin  Time 

Lat  (N) 

Lon  (E) 

1985217 

17:42:58.7 

59.3 

6.59 

1985218 

17:50:07.9 

59.3 

6.95 

1985290 

10:00:00.4 

59.3 

6.95 
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Figure  36.  Figures  Illustrating  the  noise  reduction  gain  and  beam 
loss  for  western  Norway  earthquakes  as  processed  by  beaming  and 
adaptive  filtering  with  and  without  removing  the  site  transfer 
functions. 
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FUTURE  WORK 


During  the  remaining  part  of  this  contract,  we  plan  to  apply  the  various  techniques 
we  developed  to  real  data  on  a  large  scale;  little,  if  any,  theoretical  or  code  development  is 
planned.  In  particular; 

Under  the  task  of  multi-channel  analyses  of  regional  data  we  plan  to  test  our  source 
analysis  and  azimuth  finding  approaches  described  in  this  report  on  a  large  number  of 
quarry  blast  recordings  depending  on  the  data  availability  at  NORESS,  ARCESS,  FINESA 
and  GERESS.  Since  the  acquisition  of  auxiliary  information  on  quarry  blasting  practices 
has  started  by  the  NORSAR  group,  independent  data  will  be  available  for  verifying  our 
results. 


Under  the  task  of  seismic  imaging  we  shall  apply  our  techniques  to  body  wave  data 
of  previously  studied  events,  both  earthquakes  and  explosions.  We  plan  to  apply  the 
technique  to  progressively  more  difficult  problems  starting  with  events  that  have  clear  pP 
and  sP  arrivals.  The  images  will  be  examined  for  signs  of  source  finiteness  and 
complexity.  An  automatized,  interactive  program  package  is  being  written  that  would 
enable  us  to  apply  our  imaging  technique  for  depth  estimation  and  discrimination  in  a 
manner  similar  to  the  envisaged  applications  of  Pearce's  algorithm.  This  package  can  be  a 
part  of  an  operational  discrimination  system  in  the  future. 

We  shall  resume  our  work  on  surface  wave  imaging  and  investigating  strain  release 
phenomena.  We  concluded  that  our  most  recent  efforts  were  severely  degraded  by  the  fact 
that  numerous  surface  wave  recordings  in  the  CSS  data  bases  had  erroneous  absolute 
times.  We  shall  acquire  more  data  for  which  we  have  accurate  timing  information. 
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CONCLUSIONS 


Work  during  this  repoi  ting  period  leveaied  the  great  potential  of  various  statistical¬ 
time  series  analysis  methods  for  solving  seismological  problems  relevant  to  nuclear 
monitoring. 

In  this  exploratory  study  of  regional  seismic  waveforms,  extending  the  concept  of 
spectral  factorability  of  leleseismic  body  wave  data  (Filson  and  Frasier,  1972),  we  have 
demonstrated  that  spectra  of  regional  phases  can  also  be  decomposed  into  source  and 
recording  site  factors  provided  that  the  signals  originated  from  a  limited  source  region. 
Groups  of  events  located  close  to  each  other  can  be  identified  by  the  fact  that  their 
waveforms  can  be  reconstructed  from  their  spectral  source  and  site  factors. 

In  addition,  it  was  also  found  that  events  within  such  factorable  groups  could  be 
further  subdivided,  according  to  relative  interevent  coherences  between  them.  The  events 
which  show  high  inter-event  coherence  are  probably  both  close  to  each  other  and  have  very 
similar  source  mechanisms,  albeit  different  source  time  functions.  The  opposite  must  be 
true  for  events  that  do  not  show  high  coherence.  Cross-event  equalization  filtering  between 
coherent  events  using  short  time  domain  filters  resulted  in  increases  in  waveform  similarity, 
and  no  such  increase  seemed  possible  between  events’  pairs  with  low  interevent  coherence. 

Site  transfer  functions  between  sensors  located  at  the  extreme  ends  of  NORESS 
appear  to  be  complex,  describable  only  with  transfer  functions  with  impulse  responses 
longer  than  a  second.  Assuming  the  same  degree  of  crustal  heterogeneity  for  a  source 
region  this  may  imply  changes  in  waveforms  similar  in  nature  over  small  displacements  of 
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source  location.  Further  work  needs  to  be  done  in  exploring  the  nature  of  such  transfer 
functions  for 

a)  Grouping  events  with  respect  to  relative  location; 

b)  Identifying  events  with  differing  source  mechanisms;  and 

c)  Finding  differences  in  the  source  time  functions  between  closely  spaced  events. 

We  have  derived  some  new  approaches  for  imaging  seismic  sources  and  seismic 
source  inversion.  The  methods  are  based  by  finding  the  maxima  of  statistics  that  are 
closely  related  to  those  used  in  F-K  analyses  as  functions  of  source  parameters.  The 
method  is  quite  general,  and  can  be  used  for  analyzing  both  body  and  surface  waves. 
These  algorithms  also  provide  a  means  for  evaluation  of  the  performance  ot  time  domain 
direct  modeling  and  source  inversion  methods  as  currently  practiced.  Our  analyses 
indicated  that  the  time  domain  fits  using  long  period  data  are  generally  not  precise  enough 
to  resolve  the  source  depth  to  a  few  km  and  other  spatial  details  with  the  accuracies  often 
claimed.  This  casts  doubt  on  many  results  in  the  literature.  Other  simulations  indicated  thai 
the  same  approaches  can  create  reasonably  resolved  images  of  faults  and  can  be  used  for 
the  separate  estimation  of  compressional  and  strain  release  components  of  energies  from 
explosions. 

Our  source  imaging-inversion  methodology  was  applied  to  two  deconvolved 
Kazakh  nuclear  explosions  and  the  El  Golfo  Earthquake  in  order  to  obtain  depth  estimates. 
The  method  identified  both  explosions  as  shallow,  with  depths  below  5  km,  and  the  depth 
estimate  for  the  El  Golfo  event  turned  out  to  be  14  km.  These  are  reasonable  results 
indicating  that  these  methods  can  be  applied  for  seismic  discrimination.  Imaging  Kazakh 
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explosion  sources  using  surface  waves  were  hampered  by  demonstrably  erroneous 
absolute  timing  information  in  the  CSS  data  we  retrieved. 

Applications  of  combined  adaptive  filtering  and  site-effect  compensation  processing 
schemes  to  teleseismic  (UK  array)  and  regional  (NORESS)  data  indicated  that  such 
methods  can  siginificantly  enhance  the  signal-to-noise  ratios  for  P  {Pn)  arrivals  should  such 
computer-intensive  analyses  prove  to  be  desirable  for  an  event. 
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APPENDIX  A 

SOURCE  IMAGING:  STATISTICAL  THEORY 

As  we  stated  in  the  introduction  of  this  report,  one  of  the  main  goals  of  analysis  of 
seismic  data  is  to  derive  some  properties  of  the  source.  In  nuclear  monitoring  we  desire  to 
discriminate  between  events  of  various  types,  such  as  earthquakes  and  explosions,  either 
through  the  direct  discrimination  of  source  mechanisms  or  indirectly  through  determination 
of  the  source  depth.  It  would  be  desirable  to  do  this  by  some  visual  identification  of 
sources  after  some  preprocessing.  In  conventional  processing  of  P  wave  data  such 
identification  can  be  performed  through  lining  up  the  short  period  P  waves  and  observing 
the  pP  "moveout'',  i.e.,  the  changes  in  pP  time  with  the  values  of  slowness.  Alternatively, 
one  may  apply  some  criteria  using  amplitude  measurements  in  some  time  windows 
covering  possible  times  of  depth  phases  to  ascertain  that  an  event  belongs  to  either  the 
earthquake  or  explosion  category  (Pearce,  1980).  Although  such  methods  are  clearly 
useful,  they  are  somewhat  tedious,  and  some  future  automation  of  such  approaches  is 
clearly  desirable.  In  a  broader  sense,  going  beyond  the  concept  of  point  sources,  detailed 
mapping  of  the  sources  of  seismic  radiation  is  also  very  important  in  studies  of  the 
mechanisms  of  earthquakes.  In  order  to  understand  the  faulting  processes  we  need  to 
know  how  the  fault  dislocation  developed  in  time  and  space. 

The  discussion  below  establishes  a  proposed  definition  of  the  "source  image",  i.e., 
the  estimated  power  of  seismic  sources  as  functions  of  space  coordinates  in  the  source 
region  and  derives  confidence  intervals  for  this  image.  We  consider,  as  usual,  the 
observed  time  vector  y  (t)  observed  over  N  array  elements  or  stations.  The  discrete 
Fourier  transform  of  this  can  be  written  as 

y  =fs  +  n  (1) 
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where  5  =  ifjr,  v  j  is  the  Fourier  transform  of  the  theoretical  signal  at  frequency  v  and  at 
source  position  vector  x  ,  usually  specified  as  a  combination  of  location  and  depth.  The 
vector/ =  f  (x,  V )  is  the  transform  of  the  function  that  combines  with  s  to  produce  the 
data  y.  We  can  think  off  as  some  combination  of  Green's  functions.  The  vector  /  is 
assumed  to  be  known,  whereas  the  signal  i  is  deterministic  and  unknown.  The  noise 
vector  n  is  assumed  to  be  spatially  white  (uncorrelated  among  stations)  with  power  spectral 
matrix 
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at  frequency  v ,  where  Pn(  v  J  is  the  noise  power  spectrum  of  each  channel  and  I  is  the 
identity  matrix  information  carried  by  the  observed  vector  y.  One  can  simply  plot  the 
power  in  the  generalized  transform  J*y  as  the  beampower 
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An  alternative  measure  with  higher  resolution  is  Capon's  C,  defined  as 
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which  can  also  be  written  in  the  form 


C  (x,v)=-^  1  + 
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which  shows  its  relation  to  the  generalized  beam.  It  is  convenient  to  define  the  residual 
power  as 


^  (jc,  V )  =  jy  -  B  (x,  V ) 
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so  that 


C  (X,  v)=- 
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B(x.  V) 
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the  constant  c  is  an  arbitrary  number  chosen  to  make  the  matrix  (y  y*  +  I)  non¬ 
singular. 


A  second  alternative  estimator  is  Shumway’s  F,  defined  as  the  likelihood  ratio  test 
statistic  resulting  from  testing  the  model  (1).  This  leads  to 

=  (8) 


which  is  closely  related  to  the  other  two. 

Capon's  C  and  Shumway's  F  have  similar  resolving  capabilities  since  they  are  both 
functions  of  B  and  R.  They  are  both  better  than  the  generalized  beampower  from  this  point 
of  view  since  Equation  (3)  does  not  incorporate  the  residual  noise  power  R.  We  shall 
show  below  that  the  F  is  superior  statistically  to  the  otlier  two,  bodi  for  detecting  the 
predefined  signal  and  for  providing  an  estimator  for  the  signal  image. 
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If  the  noise  is  Gaussian,  it  is  easy  to  show  that 


2  B{x,v) 
^n(v) 


(9) 


2 

where  2  denotes 


the  chi-squared  distribution  with  non-centrality  parameter 


o2  2|S(x.  V  f 


(10) 


and  ~  is  a  notation  for  "approximately  distributed  as.”  Note  that  5^  is  twice  the  "signal-to- 
noise"  ratio 


^  V )  = 


\S(x,vf 

) 


(11) 


It  is  this  signal-to-noise  ratio  that  we  propose  to  define  as  the  image.  It  is  also  easy  to 
show  that 


2Rix,v) 


(12) 


and  is  independent  of  B(x,  v).  It  follows  that 


^2,2^-2(°  1 - - 


(13) 


so  that 
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F(x,v)  =  (A/.l)f^ 


(14) 


has  the  non-central  F  distribution  with  the  non-centrality  parameter  6^  given  by  Equation 
(11).  Hence  the  distribution  of  F  depends  only  on  the  signal-to-aoise  ratio  ^{x,  v)  for  a 
given  number  of  sensors  N. 

We  may  examine  this  further  by  noting  that  for  an  F  with  ni  and  n2  degrees  of 
freedom 

=  +S^)  (15) 

'  «1.  /I2  ^2  -  2  ^ 


so  that 


£(F(x,v))=^^(l  +  |V,  V)) 


(16) 


and  the  mean  value  is  essentially  the  signal-to-noise  ratio. 


^  K  =  ?■  ,  J,  ..  1 2«i  +  45^  {n,  +5^)  ] 
(n2-2)(n2*4) 


(17) 


can  be  used  to  compute  the  variance  with 


variF)  =  E{F^)  -E^{F) 


(18) 


and  n]  =  2,  n2  =  2  N  -  2.  The  variance  depends  on  the  signal-to-noise  ratio  that  we  are 
trying  to  estimate  in  a  rather  complicated  way  and  the  normal  approximation  will  not  be 
very  good  for  obtaining  an  approximate  confidence  interval. 
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For  any  x  and  v,  however,  we  can  obtain  an  exact  confidence  internal  for 
E(F(x,v)),  8^  or  in  the  following  way.  We  know  from  (16)  that 


^2 

4 


(/V-2) 

(A/-1) 


F(x,v)-  1 


(19) 


is  an  unbiased  estimator  of  the  signal  to  noise  ratio  (11),  which  we  are  defining  as  the 
image.  LetF(x,v)  be  some  observed  sample  value  of  Equation  (14).  Let  5^;  and  5^2  ^ 

two  values  of  5^  such  that 


Pr[F. 


id,  )^F(x,v))  =  a/2 


(20) 


and 


Pr{F^^,^i&^)^Fix,v)}  =  a/2  (21) 

The  resulting  1-a  confidence  interval  for  5^  ^  ^2 ).  The  interval  for  this 

signal-to-noise  ratio  is  (1/2  ,  1/2  8^2  )•  This  interval  can  also  be  connected  into  a  (1- 

a)  confidence  interval  for  E(F(x,v)) ,  using  Equation  (16). 

The  above  procedure  can  be  used  to  get  a  (1-a)  percent  confidence  interval  for  each 
X  and  V.  Since  the  intervals  are  dependent,  for  n  of  them  one  can  only  assign  an  overall 
confidence  of  1-na  Therefore,  in  order  to  assign  small  confidence  to  the  upper  and  lower 
surfaces  of,  say,  .90,  with  100  points  on  the  surface,  each  separate  point  should  have 
a=.001  or  99.9%  confidence. 
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APPENDIX  B 


IMPLEMENTATION  OF  MULTI-CHANNEL  DECONVOLUTION 
IN  THE  X-WINDOWS  ENVIRONMENT 

In  order  to  facilitate  the  use  and  the  testing  of  the  multi-channel  deconvolution 
method,  we  have  written  a  version  that  makes  use  of  the  X-window  environment.  The 
program  is  menu-driven  with  the  following  options  for  viewing  in  the  Multi-channel 
Deconvolution  Display  (MDD): 

1)  Source  estimate  waveforms  and  resolution  kernels; 

2)  Site  ansfer  function  (time  domain)  displays; 

3)  Source  spectral  estimates; 

4)  Site  response  spectra;  and 

5)  Original  vs.  reconstructed  waveforms  with  associated  correlation  coefficients. 

Examples  of  these  various  kinds  of  displays  arc  shown  in  Figures  B1  to  B4.  For 
further  explanations  we  refer  to  the  figure  captions.  When  MDD  is  activated  it  creates  a 
window  and  waits  for  user  input.  All  input  is  entered  by  using  a  mouse.  The  right  mouse 
button  produces  a  menu  of  the  options.  The  option  is  chosen  by  moving  the  mouse  so  that 
the  desired  option  is  highlighted  and  the  left  mouse  button  is  pressed. 

Once  the  option  is  chosen,  a  window  displaying  the  desired  data  appears.  If  all  the 
data  will  not  fit  the  window,  the  first  few  waveforms  will  be  displayed  and  the  rest  can  be 
viewed  by  scrolling.  A  scrollbar  appears  at  the  left  edge  of  each  window.  It  can  be  used  to 
scroll  the  data  within  that  window  by  standard  XI 1  scrollbar  manipulations.  Each  new 
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window  will  be  created  on  top  of  those  currently  displayed.  All  windows  can  be  resized, 
moved  and  overlayed  using  standard  X  facilitites. 
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Multichannel  Deconvolution  Dlspla; 

|r«tt  butt«n  f«r  mm 
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Figure  Bl.  Display  of  source  time  function  estimates  with  the  estimated 
resolution  kernels. 


Hulticharinol  Deconvolution  Display 

rr«««  rt^  bult«n-f«r  mtu 
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Figure  B3.  Joint  display  of  source  time  function  estimates  and  their 
spectra. 


Multichannel  Doconvolutlon  Displa; 

^tit  rl#»t  button  for  m««j 


Figure  B2.  Display  of  time  domain  site  transfer  function  estimates  with 
their  spectra. 


Multichannel  Doconvolutlon  Dlspla 

right  buttwi  for  mtu 


comparing  the  data  with  their  reconstructions  from 
b1  factors.  The  value  of  the  time  domain  correlation 
each  pair  of  traces. 
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